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ABSTRACT 


Analytical solutions are obtained for steady-state symmetric and 
unsymmetric solidifications of forced laminar liquid flows in parallel- 
plate channels with convective boundary conditions. The classical 
Graetz type analysis based on the physical model of Zerkle and Sunder- 
land is made by using the confluent hypergeometric function for the 
case with uniform convective coolings at upper and lower plates and 
the case with perfect insulation at lower plate and uniform convective 
cooling at upper plate. Numerical results are obtained for liquid soli- 
dification-free length, ice layer thickness, pressure drop, bulk tempera- 
ture, heat transfer rate and local Nusselt number for a range of values 
of Biot number and superheat ratio. The configuration and thermal con- 
ditions under consideration may be encountered in heat exchangers using 
cryogenics and freezing of ice sheets on northern rivers and lakes. 

The Graetz problem with axial heat conduction effects in 
parallel-plate channels considering both the upstream and downstream 
regions is extended to the case of convective boundary condition. It 
is found that analytical solution in freezing zone cannot be obtained 


for low Peclet number flow regime and numerical solution must be sought. 
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CHAPTER I 


INTRODUCTION 


1.1 Liquid Solidification in Parallel-Plate Channels 

In recent years, the problem of liquid solidification in cir- 
cular tubes with fully developed laminar flow at the thermal entrance 
has been studied by several investigators for the thermal boundary con- 
ditions of uniform wall temperature and convective cooling. The solidi- 
fication of liquids flowing through tubes or channels is of technical 
importance and may occur in water mains, process equipment and other 
various hydraulic systems subjected to extreme ambient temperatures. 

The solidification of flowing liquid in parallel-plate chan- 
nels with convective boundary conditions is encountered in natural 
phenomena such as the freezing of northern rivers and lakes and in 
various industrial processes such as the solidifying of metal castings 
in molds, freezing of foods and solidification within liquid flow heat 
exchanges uSing a cryogen as the coolant. For the specific configura- 
tion of the parallel-plate channel, one notes that the liquid solidifi- 
cation problem has been studied so far for the case of uniform wall 
temperature condition only. 

The purpose of this investigation is to study a class of 
liquid solidification problems in parallel-plate channels with convec- 


tive boundary conditions. 


1.2 Statement of the Problem 


The problem is concerned with plane Poiseuille liquid flow in 
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parallel-plate channels with surface solidification. At the thermal 
entrance (X = 0), the liquid has a fully developed velocity profile 
and a uniform temperature distribution. In the downstream region 

(X > 0), heat is transfered by convection from the channel surface 

to the surrounding medium with an ambient temperature below the free- 
zing temperature of a warm flowing liquid. 

With a step decrease in wall temperature, ice will cover the 
entire wall starting from the thermal entrance X = 0. In constrast, 
with a uniform external convective cooling, ice will not form until a 
sufficient amount of sensible heat has been removed. Thus, the thermal 
entrance region consists of a solidification-free zone and a freezing 
zone where a growing solid layer appears along the channel wall. In 
this investigation, the liquid flow rate into the cooled section is 
assumed to be maintained constant and only the steady-state condition 
is under consideration. 

Apparently, the flow field in the freezing zone is characte- 
rized by a boundary layer flow due to the interface between the solid 
and the liquid phases moving inward. A fully developed flow is reached 
further downstream. The present investigation utilizes the physical 
model employed by Zerkle and Sunderland [1] which assumes that the axial 
velocity profile of the liquid remains parabolic even as the liquid 
accelerates due to a constriction of the flow area by the solid phase. 
This basic assumption together with the assumption that the axial vari- 
ation of the solid layer thickness is gradual leads to a considerable 
simplification in theoretical analysis and enables one to employ a 
Graetz type analysis in the freezing zone. 


For parallel-plate channels, various combinations of thermal] 
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boundary conditions at the upper and lower plates are possible. How- 
ever, the following two cases are under consideration in this study: 
1. Uniform external convective cooling at upper and lower plates 
(Chapter 2). 
2. Different but uniform external convective coolings at upper and 
lower plates (Chapter 3). 
The first case represents a symmetrical solidification problem and 
the second case deals with an unsymmetrical one. For the second case, 
numerical results are obtained only for a specific case with perfect 
insulation at the lower plate and a uniform convective cooling at the 
upper plate. The analytical solutions for the foregoing two cases are 
obained by the eigenfunction expansion method using the confluent hy- 
pergeometric function. 

For low Peclet number flow regime, axial heat conduction 
effects are known to be important in thermal entrance region problem. 
Thus, it is of practical interest to consider the axial conduction 
effects on liquid solidification in parallel-plate channels with con- 
vective boundary conditions (Chapter 4). It is found that analytical 
solution in the freezing zone cannot be obtained when one considers the 
axial heat conduction effects. Consequently, numerical results are 
obtained only for the solidification-free zone and are presented in 


Chapter 4. 


1.3 Format of Presentation 
A remark regarding the structure and the method of presen- 
tation for the present thesis is now in order. In view of the analytical 


solutions obtained for three different cases, it is decided that each 


tybug2 ‘2tds nt no J 


‘eenula gnc bom ipa i 


5,4 
ee ‘ 


bins Yeqgu 36 a gvizoevnoo ee. 

. (8 | 

oes setdow vclseattibtfoe.. fsotvdemne. a0 
4na> beeoer an! ~od Gene Paatadommeeny 6 ing 
dastraa fifty 9262 200foedt & 107 eine ben tside ‘aie 2d fuze .y ai 
ott 46 portveas avia 28 ¥naD eroticu # bos stata Newol aid 36 6 tof 
oii 29@n3 owt oni opetot Sfi iw), enoltiuioe Fiigeten, oT ten 
“vi Inaw!tnap ot? antay bottem ablapeqns get) oautnapts ans a hentat 
moFtonwt sindemsaeg” 

notgovbns geet fetes .cateey wolt tadmin tufoad wal 964 ; 
mefotri| eng as rete \owrert! ot jJasiveqet od af Wand ave eiaeite 
“ol topeeeo. fetse sé) wabiencr at tzssetnt Isotioesd We af a .2udT 

- i s{heruts epeta-feffoeq of soresattibifoe eipehl ne eaoette. 
loatsy(eme god! HoeRaT gL .1¢ ve3 dad) senolatbhos qabnwad av itoey a 
aft! easpl race enc cate Benleies ad Jornis9 anas oni aes? era ob notduloe 
et atfvet> [sate <Tineugeenc) .2jootts norsoybn0s +404 eta 

ni bolnstend. sa Gms gues cort-nolsaati tht foe eng yeh Ulan bontesdo 
ee Poop 


notisinszsss 2038 

-nees6g Io bwts od Re ts arutouy2 sd? pabbrspe7 seme A 
[sotivlens off Ye woly of . Sein? won et atenaa, anes 1] 86 u 
fee. fale hoblaeb etd Aaah sneneTAtnigantdy bent: ad ‘ 


Te % ; 2 as 


| 


ad ay, f 
: .) 


wee PS 


case will be treated independently in each Chapter. Thus, the litera- 
ture review can be made separately for each case and the scope of re- 
sults, conclusions and significance for each case can be examined 


separately. 
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CHAPTER I] 


LIQUID SOLIDIFICATION IN A 
CONVECTIVELY-COOLED PARALLEL-PLATE CHANNEL 


A theoretical analysis of liquid solidification in the thermal 
entrance region of a parallel-plate channel with a uniform external 
convection is carried out by using the confluent hypergeometric func- 
tion in the solution for laminar flow and steady-state conditions. 
Theoretical results are obtained for liquid solidification-free length, 
liquid-solid interface profile, pressure drop, bulk temperature, heat 
transfer rate and local Nusselt number for the cases of Biot numbers 
Bie = e250. cle ee Superneat Yatios.c = 1, 2. 4, 8. Pr= 13,2 
and Kp/k. = 0.25 using 20 eigenvalues. The results can be used in 
assessing the effect of thermal insulation and predicting the thermal 
insulation required to protect a channel flow from a severe cold en- 


vironment. 
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NOMENCLATURE 


parameters in Kummer's equation 

Biot numbers, Not /k and (ky /k, )Bi 

constants in Eq. (8) 

specific heat at constant pressure 

eigenconstants, Eqs, (3) and (29) 

functionvon xs"Eqv'(25) 

local heat transfer coefficient, Eqs. (55) and (56) 
overall heat transfer coefficient between inner wall and 


ambient defined by -k(9T/9Y) Tah or =U 


Y=L_ “0''y= 
thermal conductivity of liquid 
half-distance between parallel plates 
confluent hypergeometric function 
eigenfunctions, Eqs. (26) and (4) 
Nusselt number 

pressure and pressure at X=0 

Peclet number, 4 U t/a 

Prandtl number, v/a 

dimensionless pressure drop, 2(P-P)/pUS 


total heat transfer rate and dimensionless Q, Eq. (50) 


local heat transfer rate, Eqs. (55) and (56) 


= liquid temperature, bulk temperature and freezing tem- 


perature. 
uniform entrance temperature and ambient temperature 
axial and transverse velocities 


average axial velocity 
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Greek Letters 


Subscripts 


aS 


Teez 


solution of Kummer's equation 
rectangular coordinates 
axial coordinate with origin at solid-phase entrance 


(X/LPe), (Y/L) 


= (X'/LPe), (Y/L) in freezing zone 


dimensionless solidification-free length 


transformed variable 


thermal diffusivity 

eigenvalues in Eqs. (3) and (29) 

transverse coordinate of liquid-solid interface and 6/L 
superheat ratio, (Ty - Te)/(T, - T.,) 

dimensionless transverse coordinate, y/6 = Y/é 
dimensionless temperature difference, (T - T,,)/(T9 = aale) 
dimensionless bulk temperature, Eqs. (47) and (48) 
CLG oe } 

dynamic viscosity and kinematic viscosity 

density 


dimensionless temperature difference, (T - T,)/(Tg - Te) 


liquid and solid-phase 


solidification-free and freezing zone 
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coe eintroduction 

The problem of liquid solidification and pressure drop in a 
circular tube with laminar flow has been studied by Zerkle and Sunder- 
land [1], DesRuisseaux and Zerkle [2], Bilenas and Jiji [3,4] and Hwang 
and Sheu [5] for the case of a uniform wall temperature below the free- 
Zing point. When a pipe is exposed to a meteorological environment or 
is buried underground in the permafrost region with thermal insulation, 
then the convective boundary condition (thermal boundary condition of 
a third kind) outside the pipe is more realistic. The case of convec- 
tive boundary condition was studied by Zerkle [6] and Lock, Freeborn 
and Nyren [7]. One notes that the uniform wall temperature represents 
a limiting case of the convective boundary condition with Biot number 
equal to infinity. 

The solidification problem of liquids flowing through ducts 
or channels is of considerable practical importance in various engineer- 
ing processes or devices ranging from casting of metals in metallurgical 
processes to freezing of water pipes under northern climatic conditions. 
Lee and Zerkle [8] presented theoretical results for liquid solidifica- 
tion in a parallel-plate channel with laminar flow using the physical 
model of Zerkle and Sunderland [1]. 

The present problem is concerned with the solidification and 
pressure drop of the superheated laminar liquid flow in a parallel- 
plate channel subjected to a uniform convective cooling condition with 
the surrounding sink temperature below the freezing point. The freez- 
ing in parallel-plate channels occurs, for example, in various plate 
coolers. The analysis assumes steady-state conditions [1,2,6] and a 


constant finite thermal resistance between the channel inner wall and 
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the external cold environment. With a convective boundary condition, 
liquid soldification will not occur over an initial thermal entrance 
length due to insufficient cooling. Thus, the thermal entrance region 
consists of a solidification-free zone and a zone with solid-phase 
where soldification grows inwards with the distance along the channel. 
The purpose of this investigation is to study the effects of superheat- 
ing and thermal insulation on heat transfer, solidification-free zone, 
liquid-solid interface profile and pressure drop in parallel-plate chan- 
nels. The solution for thermal entry problem is obtained by the eigen- 
function expansion method using the confluent hypergeometric function 


[9-12]. 


2.2 Problem Formulation and Physical Model 

Consider the steady-state freezing of a liquid with constant 
physical properties flowing between two parallel plates from the up- 
stream adiabatic section to the downstream section with a uniform con- 
vective cooling (see Fig. 1). It is assumed that the flow is a steady 
plane Poiseuille flow with a uniform entrance temperature Tg at X = 0 
which is greater than the freezing temperature Te. The basic assumption 
[1,6,8] that the axial velocity profile of the liquid in the region with 
surface solidification remains parabolic even as the liquid accelerates 
due to a continuously increasing constriction of the flow area along the 
channel axis proves to be a good approximation [3,5] and leads to a 
considerable simplification enabling one to employ a Graetz type analy- 
Sis. This assumption together with the assumptions of negligible vis- 
cous energy dissipation and free convection will be used in the present 


analysis. An order of magnitude analysis given in appendix 1 reveals 
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that axial heat conduction may be neglected if 1/Pe << 1. 


2.3 The Liquid Solidification-Free Zone 
2.3.1 Analysis: 
The energy equation and its boundary conditions describing 


the temperature field in the liquid are, 


oT OT 
U— = Q——— 

OX yye 
T(0,Y) = Tp 


91(X,0)/aY = 0 
-kaT(X4L)/3¥ = hol T(X.L) - TJ 


The Poiseuille velocity profile which satisfies the constant mass flow 


rate criteria fe Udy = LU. is given by 


3 Y\2 
u=3u 0 - (4 


Thus, the energy equation in dimensionless form and the ther- 


mal boundary condition can be written as 


2 
3 2,00 0.0 
a es ae fd) 
8 OX ay" 
6(0,y) = 1, 90(x,0)/ay = 0, 96(x,1)/ay = -Bid(x,1) (2) 


where x = X/LPe, y = Y/L, 9 = ae saline le) and Bi = hl /k and all 
other symbols are defined in Nomenclature. By separation of vari- 


ables, the solution to Eq-(1) is found as 


9 = jE 1E5Y; (yexp(- 30x) (3) 
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where the eigenvalues ot, and eigenfunction ve Satisfy the following 


Sturm-Liouville system of equations: 


d“y. Hae 
dy” % (1-y ae =y0 (4) 
dY (0)/dy = 0, d¥ (1 )/dy = ~BiY <1) (5) 


In order to overcome the difficulty in obtaining the higher eigen- 
values, the series solution is replaced by a solution involving the 
confluent hypergeometric function [9-12]. By introducing the trans- 


formations, z = aay? and w(z) = exp(z/2) > ¥ 5(y). Eq. (5) becomes 


which is seen to be Kummer's equation, 


dew 
aay it (b - zs - aw = 0 (7) 
dz 


with a = aA. b = 1/2. The general solution of Eq. (7) can be 
expressed in terms of the confluent hypergeometric function [11] 


M(a,b,z) and the eigenfunction Y 5(y) becomes 


4 2 2 
¥ 5(y) = exp(-o 5 /2)[B, Mla; »1/2,0.y ) 
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where By; and Bo; are the constants. 
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2.3.2 Eigenvalues and Eigenfunctions 


From Eq. (8), one obtains 
dY.(y)/dy = exp(-o..y°/2)[-B ocyM(a.s1/2,0<y°) 
J J ta PSs J J 


+ 8B pin Me é 
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+ (4/3)B ey? (a +1/2)M(a,+3/2,5/2.0,Y°)] (9) 


va 


Applying the symmetry condition dy :(0)/dy = 0, one obtains B 0 


Opi 


indicating that all odd terms vanish. With B,. = 0, the constants 


2J 


B can be absorbed by the eigenconstants E. INeEGRE an LbUS sees. 


1j 
(8) and (9) become 


Y;(y) = exp(-a,¥°/2)N(aj41/2.0,¥°) (10) 
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The confluent hypergeometric function M(a,b,z) is 


a(at+1)...(atn-1) Ze 


Gearlog b(bt1) 2! °° °° * b(be1)-++(ben-1) ml 


With b = 1/2 and a. (1-o,)/4; one obtains 


(l-a.)+++(4n-3-a.. ) 
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Substituting Eqs. (10) and (11) into the convective boundary condition 


in Eq. (5), one obtains 
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after applying the following recurrence formula. 
zM(a,b+1,z) = bM(a,b,z) - bM(a-1,b,z) (13) 


It is convenient to express Eq. (12) in series form as 


n : 
3 + 2 = ¢ = 
O (Bi n a) 0 (14) 


The eigenvalues a, are the roots of Eq. (14) and can be found by tabu- 
lating the value of the left-hand side of Eq. (14) versus a sequence of 
closely spaced values of ae One may thus inspect the whole spectrum 
of the eigenvalues before they are improved by the variable secant 


method. 
C.o8o Eigenconstants 


Utilizing the orthogonality relationship, the eigenconstants 


E. can be determined from 
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The first twenty eigenvalues, derivatives and constants are computed 


for n = 150 and are given in Table 1 for Bi = 0.25, 0.5, 1, 2, 10. 


2.3.4 Liquid Solidification-Free Length 

The solidification starts at the point where the wall tem- 
perature reaches the freezing temperature. Using Eq. (3), the liquid 
solidification-free length X~ Can be evaluated from 


20 
2 a 2 


where Fe = (T. - T,,)/(T, - ie) For the determination of Xe iter 


convenient to define a superheat ratio € as 


en eM) (ied) 8 ek or (18) 


It is readily seen that Xe = fF(eeBia: 


2.4 The Freezing Zone 
2.4.1 Analysis 
Assuming that the axial velocity profile remains parabolic 


[1,8] in the freezing zone (see Fig. 1 for coordinate system), the ve- 
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locity components which satisfy continuity relations and boundary con- 


dition can be shown to be 


U(X'LY) = (3LU_/28)[1 - (¥/8)7] | (19) 


WGXae 


(3/2) (LU, ¥/8“)(d8/aX")[1 - (¥/8)2] 
Substituting Eq. (19) into energy equation, U(aT/ax') + V(aT/aY) 


= o(9°T/3Y°), one obtains 
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Defining the dimensionless variables, 


<I 
i} 


(Xi 7a eyeas Vee On= 6/L, o = (T-Te)/(Ty-Te)» Eq. (20) becomes 
“ ae 2 
Seok lee. ()6T [Sere SS Abe 25 (21) 
6 6 OX 6 dx Oy oy 


The boundary conditions are 


6(x,5) = 0 (interface), 36(x,0)/ay = 0 (symmetry), 


(22) 
o(0,y) = [(1 + e)/elo(x,sy) - (1/e) (entrance condition) 


where the temperature distribution at the solid-phase entrance 8(xpsy) 
can be determined from Eq. (17). Introducing the variable transforma- 
tion [1,8] from 6(x,y) to 6(xsn) with n = y/é = Y/6, one obtains the 

following system of equations which is of the same form as the classi- 


cal Graetz problem. 
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Si linahol walled, eC 23 
ax 6 (1-n°) ané Ce 
(x51) = 0, 36(x,0)/8n = 0, o(0,n) = [(1te)/e]e(xe,y)-(1/e) (24) 


Separating the variables by putting » = G(x)N(n), one obtains: 


BeeydGee 2s ed te u 
=) fest BG = 0, G(0) = = 8(xpsn) - 5 (25) 
xX 
2 
d N ° ? fe 
—zl + 62(1 - n°)N, = 0, dN.(0)/dn = 0, N.(1) = 0 (26) 
a j j j j 


Integrating Eq. (25), one obtains 


6(%) = exp [-(8/3)85/5 (27) 


Since Eq. (26) can be transformed into Kummer's equation, the solution 
is found to be 


8 £ 
t 2 1 A al 2 
N5(n) = exp(-Bn /2)M(q I 1/258. ) (28) 


Thus, the series solution for o(x,n) is 


o(on) = JF, KNy(ndexol(-8/3)85/9 & (29) 


Where K; are the eigenconstants. 


2.4.2 Eigenvalues and Eigenconstants 


The eigenvalues Bs can be determined by using the boundary 
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This leads to 
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NG ~ gle 25 8,) 


= 


(30) 
The roots representing the eigenvalues can be determined by the proce- 
dure uSing variable secant method discussed earlier 

Applying the entrance condition for o( 


sn) and using the 
orthogonality relationship, the eigenconstants can be evaluated from 


K = ihe (1-n? JN (n)9(Osn)dn/f, (1-n°) “(n)dn (31) 
20 
where 6(0,n) = wey 


The expressions for the derivatives ON (1)/93B 
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It can be shown (see Appendix 2) that the expression for K. reduces to 
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The first twenty eigenvalues and the corresponding eigenconstants 


together with dN, (1)/dn % dNog(1)/dn are listed in Table 2. 


2.4.3 Solid-Phase Layer Profile 
Considering an energy balance for steady state at the 
ligquid-solid interface, one obtains the following interface equation 


which is a characteristic feature of the solidification problem. 
Kien, (X=,0)70¥" = k aT (X Po) ak (35) 
URES ERIC I) TO ae ey Te)/(L8)]96(x51)/an (36) 
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The temperature in the solid-phase layer is described by the 


following system of equations. 


Gal dT _(X',L) 
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The temperature solution is 


Bi (T,-T,,) 
ee ee ear ee 


where Bi. = hol /k, = (k,/k, )Bi. 
The linear temperature gradient for solid layer is 
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From Eq. (35), one then obtains 


ag (14Bi .)9¢(x,1)/an 

6(x) = Bi [80(%,1)/on-(k 7k e)] (40) 
E 20 dN.(1) Sink 

where dx,1)/9n = pak if expL(-8/3)85 Jo = oe (41) 


Solutions for 6 (x) can be found numerically by a method of successive 


substitution[8]. 


Mae Pressure Drop and Heat Transfer Characteristics 
2.5.1 Pressure Distribution 

| The pressure drops linearly with X in the solidification- 
free zone with a pressure gradient dP/dX = (-3uU_)/L°. Introducing 
the dimensionless pressure difference p = (Py - P)/(pU“/2) and upon 


integrating, one obtains 
p = 24 Prx, —Ox<x¢ (42) 


To determine the pressure drop in the freezing zone, one 
integrates the axial momentum equation, UdU/dX' + VOU/dY = -dP/pdX' 


+ \ (a¢usax'® + 3¢u/ay’), from Y = 0 to Y = 6 and obtains 


@ 
Guercne 6 dP § ay U8 
axe Jo Y= - 5 axe * Wo Syeda Wee ay (43) 
: : : ave 2 a2 
Substituting U from Eq. (19) and assuming [1] (dé/dX')”, (d 6/dX'")<<1, 


Eq. (43) in dimensionless form becomes 
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dx Sse dx x3 
After integrating from x = 0 to x = x and noting that p = 24 Pr Xp at 
xX = 0, one obtains 

vie ee ALR x dx 
ACG) 5 (2 1) + 24 Pr(x¢ + Jo x (45) 


One notes that the first term represents the momentum effect and the 


second term is due to the viscous effect. 


2.5.2 Bulk Temperature 
The dimensionless bulk temperature in the solidification- 


free zone is 
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Using Eq. (4), one obtains 
rah be 
1(4_yé mi poel eae 
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Thus, Eq. (46) reduces to 
e, = "= > 52] E = expl(-8/3)a ix J/a, O<x<X ¢ ) 


Following the same procedure, the dimensionless bulk temperature in 


the freezing zone can be shown to be 
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After making use of Eq. (26), one obtains 
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255-5) neat lranster 
The dimensionless total rate of heat transfer from the 
liquid per unit depth of the channel in the solidification-free zone 


from x = 0 to x = x can be found from 


Q, = Q/2LU_oc, (T)-Te) (50) 


p 


where Q = 2f-kdT(X,L)/3¥*dX, — OSX<xe. 


One notes that 0 assumes a value of 1 when the liquid is cooled to 7ts 


freezing temperature. Using Eq. (3), one obtains 
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Similarly, the dimensionless heat transfer rate in the freezing zone 


from x = 0 to x = x is found as 


dx (52) 


Thus, the total heat transfer rate between the thermal entrance (x=0) 


and any axial location (x=x) in the freezing zone is simply given by 
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2.5.4 Local Nusselt Number 


The local Nusselt number provides further insight in under- 


Standing the heat transfer mechanism for the present liquid solidifica- - 


tion problem and is defined by 


Nu = h(2L) = =f 36 (x51) /, ; 0< x< x (53) 
ky ay b f 


Nu = = - EBON Tg Bey] Xp (54) 
6 


In deriving the above equations, one notes that 


q=- k, Te} = SU Fe O<X<X ¢ (55) 
SL) s X>X ¢ (56) 


2.6 Results and Discussion 

The theoretical results for the solidification-free length 
X~ Versus Superheat ratio € are shown graphically in Fig. 2 with Bi 
= 0.25, 0.5, 1, 2, 10. As expected, for a given value of e, Xe in- 
creases with the decrease of Biot number. The limiting case Bi =o 
corresponds to a uniform wall temperature with Xe = O. The superheat 
ratio € = 0 v 10 is considered to be the practical range in engineering 
applications [7]. 

Plots of the solid-phase layer thickness 5 versus the en- 


trance distance x are given in Fig. 3 for € = 1, 2, 4, 8 and Bi = 1, 
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2, 10, 100, ~. The location where the departure of the curve from 6 


= 1 occurs represents the solidification-free length x The solid- 


¢ 
phase increases its thickness with x but stops short of reaching 6 
= 0 since a constant mass flow rate is assumed in the analysis. This 
precludes the "freezing-shut" of the channel. The effects of Biot 
number or superheat ratio on the solid-phase layer profile can be 
seen Clearly in Fig. 3. The limiting case Bi = ~ corresponds to a 
constant wall temperature studied by Lee and Zerkle [8]. An examina- 
tion of Eq. (40) reveals that as Bi > ~, the expression given by Lee 
and Zerkle [8] is recovered. 

The pressure drop results are presented in Fig. 4 for Pr 
= 13.2 (water). The straight line part of the pressure distribution 
curve represents a pressure drop in the solidification-free zone which 
is independent of Bi and ec. When ice is being formed, the flow acce- 
lerates due to the reduction of the flow area and the pressure drop 
increases rapidly in the freezing zone. Thus, a large pressure dif- 
ference is required to maintain a constant flow rate when the ice 
growth nearly fills up the channel. Eqs.(42) and (45) show that the 
pressure drop increases with Prandtl number. 

The dimensionless bulk temperature distributions for Bi 
= 1, 2, 10,” are shown in Fig. 5 with ce = 1 and 4. Noting that 8 
= 1/(lte) in Eq. (48), one sees that the asymptotic value for 6. 


is 6,. When 6 the bulk temperature approaches the freezing 


7 b °F 
temperature. For a given value of e, the thermal entrance length 


for approaching an = O decreases with the increase of Biot number. 


The amount of sensible heat to be removed before the bulk 


temperature approaches the freezing temperature is of considerable 
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practical interest. The dimensionless heat transfer rate 0 as a func- 
TLOn On xis) Shown MN Pigewoeton bis=s, 2.4) 10, © with o=- lvand-4; 
It is noted that 6 > 1 as 6, > (see Fig. 5) and the effect of super- 
heat ratio is considerable. 

The local Nusselt number results are shown in Fig. 7 and the 


results may be contrasted to the axial variations of 6, and 0. It is 


b 
of interest to note that a local minimum value for Nu exists. The 
effect of the Biot number and superheat ratio is evident. In this study, 
the solution is terminated when 6 < 10°? is reached since the physical 
model may not be valid with further reduction of the flow area. One may 
wish to define a local Nusselt number in the freezing zone by using the 
freezing temperature Te instead of the surrounding fluid temperature ie 
in Eq. (56). Then the discontinuity in local Nusselt number occurs at 
the end of the solidification-free zone and the results are presented 

in Fig. 8. Figs. 7 and 8 show that the local Nusselt number decreases 
and increases monotonically in the solidification-free and freezing 
zones, respectively. The behavior of Nu in the solidification-free 

zone is well understood and is due to the entrance effect caused by 

the growing thermal boundary layer. For Bi = 1 and 2, the end of the 
solidification-free zone corresponds to the fully-developed condition. 
The increase in Nusselt number in the freezing zone is caused by the 
increased axial velocity due to the presence of the increasing solid 
phase layer. Further insight regarding the local Nusselt number be- 
havior can be gained by noting the following relationship and con- 


sidering the relative magnitudes of the numerator and denominator. 
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The local Nusselt number behavior for the freezing problem in a circu- 
lar tube or channel does not seem to have been elucidated in earlier 
Studies. The increase of the pressure drop in the freezing zone (see 


Fig. 4) is clearly related to the corresponding increase in local 


Nusselt number. 


ey, Concluding Remarks 

The solution of the steady-state liquid solidification pro- 
blem in a parallel-plate channel flow with uniform external convection 
is obtained by using the confluent hypergeometric function. The first 
Six eigenvalues in the solidification-free zone for Bi = 1, 2 agree 
with the available results [12] in the literature. The use of the con- 
fluent hypergeometric function in finding the eigenvalues from Eqs. 
(14) and (30) proves to be very efficient computationally and requires 
considerably less computing time than the Runge-Kutta method, for 
example. 

Recent theoretical and experimental results of Hwang and 
Sheu [5] for the liquid solidification in combined hydrodynamic and 
thermal entrance region of a circular tube with a uniform wall tem- 
perature show that the physical model used by Zerkle and Sunderland 
[1] is a reasonable one. Thus, one may reason that the present ana- 
lysis for liquid solidification in a parallel-plate channel using 
the same physical model represents a good approximate solution. The 


limiting case of the present analysis with Bi = ~ corresponds to an 
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isothermal wall studied by Lee and Zerkle [8]. 

For liquid flow between horizontal parallel plates, the free 
convection effects will appear only after the onset of a secondary flow 
in the form of longitudinal vortex rools [13,14]. The thermal instabi- 
lity problem in the solidification-free zone represents a separate pro- 
blem. However, the free convection effects in the freezing zone are 
not expected to be important since the temperature difference across 
the potentially unstable layer between the midplane and the upper liquid- 
solid interface is very small. This is in contrast to the case of tube 
flow where the free convection effects always [1] exist and are repre- 
sented by Grashof number. 

The pressure drop results shown in Fig. 4 together with the 
other results are useful in estimating the conditions [2] under which 
a parallel-plate channel may freeze shut by noting that the solution 
is terminated in this study when 6 = 10° is reached. One notes that 
solution using the confluent hypergeometric function can also be applied 
to the internal axisymmetric free-boundary problem involving freezing 
in a cooled circular tube [1-7]. More eigenvalues and eigenfunctions 
are required to obtain an accurate solution near the thermal entrance 
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Eigenvalues, Eigenconstants and Related Derivatives for 
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Table 2. 


Eigenvalues, Eigenconstants and Related Derivatives 


for Freezing Zone with « = 1, 4, Kp /k, = 0.25 
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CHAPTER III 


ASYMMETRIC SOLIDIFICATION OF FLOWING LIQUID IN 
A CONVECTIVELY COOLED PARALLEL-PLATE CHANNEL 


In this Chapter, an analysis is made of steady-state unsymmetric 
liquid solidification of forced laminar flow in a parallel-plate 
channel with different uniform external convection coolings at the 
upper and lower plates. The classical Graetz type analysis is made 
by using the confluent hypergeometric function. The case of one 
plate with perfect insulation and the other plate with a uniform 
external convection cooling is studied in detail. Numerical results 
are obtained for liquid solidification-free length, ice layer thick- 
ness, pressure drop, bulk temperature, heat transfer rate and local 
Nusselt number for Bi, Sau, Bi, ecw l00,.ee= Oe 10S Pr = 
13.2 and kK /k. = 0.25 using 10 eigenvalues. The configuration and 
thermal conditions may be encountered, for example, in heat exchan- 
pire using cryogenics and freezing of ice sheets on northern rivers 


and lakes. 
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NOMENCLATURE 


parameters in Kummer's equation 

Biot numbers at lower and upper plates, hjL/k; hoL/k 
(ko/k)Bi, 

constants in eq. (8) 

constants 

Specific heat at constant pressure 

eigenconstants, eqs. (4) and (33) 

function of x only, equation (31) 

local heat transfer coefficient, eas. (54) and (55) 
overall heat transfer coefficients defined by k3T(X,0)/aY 
= h,LT(X,0)-T, J, - kdT(X,L)/3Y = hoLT(X,L)-T] 
thermal conductivity of liquid 

distance between parallel plates 

confluent hypergeometric function 

eigenfunctions 

Nusselt number, hL/k 

Pressures at X and X = 0 

Peclet number, 4 U L/a 

Prandtl number, v/a 

dimensionless pressure drop, (P.-P)/(oU</2) 

total heat transfer rate and dimensionless Q, eq. (51) 
local heat transfer rate, eqs. (54) and (55) 

liquid temperature, bulk temperature and freezing 
temperature 


uniform entrance temperature and ambient temperature 
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axial and transverse velocities 


average axial velocity 

solution of Kummer's equation 

rectangular coordinates 

axial coordinate with origin at solid-phase entrance 
(X/3LPe), (Y/L) 

(X'/3LPe), (Y/L) in freezing zone 

dimensionless solidification-free length 

transformed variable 

thermal diffusivity 

eigenvalues 

transverse coordinate of liquid-solid interface and 6/L 
superheat ratio, (T,-T,)/(T,-T,,) 

dimensionless transverse coordinate, y/6 = Y/é 
dimensionless temperature difference, aU, ueinade) 
dimensionless bulk temperatures, eqs. (48) and (50) 
(TT I(T gM) 

dynamic viscosity and kinematic viscosity 

density 


dimensionless temperature difference, (iat) / (lead) 


liquid and solid-phase. 
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3.1 Introduction 

The freezing of liquids in forced laminar flow inside 
circular tubes is of technical importance and has been studied by 
Zerkle and Sunderland [1], DesRuisseaux and Zerkle [2], Ozisik and 
Mulligan [3], Bilenas and Jiji [4, 5] and Hwang and Sheu [6] for 
uniform wall temperature condition and by Zerkle [7] and Lock, 
Freeborn and Nyren [8] for uniform external convection cooling. The 
analysis of Zerkle and Sunderland [1] employs the basic assumption 
of a parabolic axial velocity profile even in the presence of steady- 
State ice layer and reduces the problem to Graetz type analysis which 
proves to be a good approximation confirmed by numerical solutions 
[4, 6] and experiments [1, 6]. Lee and Zerkle [9] extended the analy- 
Sis to steady-state liquid solidification in a parallel-plate channel 
with uniform wall temperature below the freezing temperature. 

Turning to the one-plate problem, the transient freezing of 
a liquid flowing over a flat plate was treated by Lapadula and Mueller 
[10] and Beauboeuf and Chapman [11] for the case of constant wall 
temperature and by Miller and Jiji [12] for the case with constant 
heat removal. Savino and Siegel [13, 14] carried out experimental 
and analytical investigations on transient solidification of a warm 
liquid flowing over a convectively cooled flat plate. Their inves- 
tigations where motivated by the application in liquid-liquid heat 
exchangers employing the coolant at a temperature below the freezing 
point of the warm liquid [13]. One notes that their test apparatus 
corresponds to the configuration with liquid flow between two par- 
allel plates. Since the convective heat transfer decreases with 


length along the plate, the ice thickness increases in the downstream 
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direction and the ice-layer profile (see Fig. 6(b) of [13]) resembles 
that of Zerkle and Sunderland [1]. 

The purpose of this study is to analyze the steady-state 
two-dimensional asymmetric solidification of laminar forced liquid 
flow in a parallel-plate channel with different, uniform external con- 
vection coolings at the upper and lower plates. The investigation was 
motivated by the applications in solidification of metal castings in 
molds, freezing of northern rivers with ice cover, and solidification 
within liquid flow heat exchangers using a cryogen as the coolant 
such as plate coolers [13]. 

One notes that the two-plate problem represents a generali- 


zation of the one-plate problem since flow over a flat plate can be 


recovered when the distance between the two plates approaches infinity. 


On the other hand, the present investigation can be regarded to be 
the extension and generalization of the analysis by Lee and Zerkle 
[9] to the convective boundary conditions at the upper “nd lower 
plates. For the convective boundary condition, the limiting cases 

of Bi = 0 and ~ represent the perfect insulation and the uniform wall 
temperature, respectively. The case Bi # O represents a finite ther- 
mal resistance between the inner wall and the cold ambient. The 
present analysis is valid for any combination of Biot numbers at 

the upper and lower plates but numerical results are obtained for 

Bi, = 0 (perfect insulation) at lower plate and Bi, = ee 10. 00 
at upper plate. This type of asymmetric solidification problem arises 
in freezing of northern rivers, discharge of warm water from power 
plants to the frozen northern lakes and various plate coolers. 


For Graetz problem in parallel-plate channels with unsymme- 
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tric boundary conditions, one may mention the works by Schenk [15], 
Dennis and Poots [16], and Cess and Shaffer [17, 18] for reference. 
In this study, the solution to the thermal entry problem is obtained 
by the eigenfunction expansion method using the confluent hypergeo- 


metric function [19, 20]. 


Bac The Liquid Solidification-Free Zone 
3.2.1 Physical Problem and Analysis 

The physical model is depicted in Fig. 1. The warm liquid 
with an established plane Poiseuille flow enters the thermal entrance 
(X = 0) at a uniform temperature TO: Different, but constant overall 
heat transfer coefficients, hy and hos between the inner wall and the 
cold ambient are prescribed at the lower and upper plates, respectively, 
in the cooling section. The liquid is cooled as it flows through the 
channel and eventually, at some downstream location, the freezing tem- 
perature is reached at the wall and solidification begins. With h, # 
hos the location of point of incipient solidification is different at 
the lower and upper plates. Thus, the thermal entrance region consists 
of a solidification-free zone and a freezing zone wnere the solid layer 
increases its thickness monotonically along the channel wall. It is 
assumed that the physical properties are constant and axial heat con- 
duction, viscous energy dissipation and free convection are negligible. 
The fully developed velocity profile is given as 
U(X.¥) = 6U,,(Y/L)[1-(Y/L)] (1) 
The energy equation in dimensionless form for the steady state and the 


thermal boundary conditions can be written as 
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Z 
99 _ 9°6 
LAU EAI ees ary | (2) 
dy 
6(O,y) = 1, 36(x,0)/ay = Bi,6(x,0), 36(x,1)/ay = -Bi,6(x,1) (3) 
her er mee (oiPe) say av /eence e(- Ty / (1-1) Bi, = hy L/k 
Bi, = hoL/k 


and all other symbols are defined in Nomenclature. By using the method 


of separation of variables, the solution takes the form 


6(x,y) = _ EY, (y)exp(-a5x) (4) 


where & and ue are the eigenvalues and eigenfunctions which satisfy 


the following system of equations. 


dy. , dY.(0) dy .(1) 
Meena ie ol Sr peey (ny) le pave] 5 
ae aL <y( y) j ag i, ) ay i, ! ) (5) 


3.2.2 Eigenvalues 
In order to overcome the difficulty in obtaining the higher 
eigenvalues, the series solution is replaced by a solution using the 


confluent hypergeometric function [19, 20]. By introducing the trans- 


formations, z = a,(y-1/2)" and viz) = ¥ sy) - exp(z/2), equation (5) 
becomes 

acy dv 04 

z—z + (1/2-z)=— - (1/4 - =e) v= 0 (6) 
dz° dz 16 


which is recognized to be Kummer's equation, 


Zand 79 br av = 0 (7) 
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with a = (1/4)-(a,/16), b = 1/2. Noting that the general solution 
of equation (7) can be expressed in terms of the confluent hyper- 
geometric function M(a,b,z) [20], the expression for Y 5(y) takes 


the form 


Y jy) = expl-a,(y-1/2)°72]: [C, sa, 1/2,0,,(y-1/2)*) 


+ Gee bout 1/2)M(a,+1/2,3/2,0.,(y-1/2) a9 (8) 
where ak and Co; are the constants. The series expansion for 


M(a,b,z) is [20] 


a(at])***(atn-1) z 


M(a,b,z) = 1 + b(btl) 21 b(btl)***(bén-1) nt (9) 


Thus, one obtains 


Y.(y) = exp[- “a. (y-1/2) £72] 
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(3-a./4)***(4n-1-0./4) 
1/2 @ j j 
+ a! nz] (ont)! al (y-1/2)"" 


-[-a,,(y-1/2)"42n41]3 (11) 


where Co is being absorbed by the eigenconstant in equation (4) and 
C. is the new constant. Substituting equations (10) and (11) into 


the convective boundary conditions in equation (3), one obtains 
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By adding and subtracting the above two equations, one can solve for 
C. and obtain the equation for O15 The subtraction of equation (12) 


from (13) yields 
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while the addition of the two equations gives 
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Equating equations (14) and (15), one obtains 


(3-a./4)**+*(4n-1-a./4) 
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(1-a0./4)+***(4n-3-0./4) a. 


[1+ by (2n)! Gy] = 0 
(16) 
In present study, only the specific case of Bi, = 0 (perfect 
insulation at bottom plate) will be studied in detail. With Bi, = 0, 
equation (16) becomes 
9 wo | (37/4) ***(4n-1-a./4) oe 
COD ih OB ana ie coecemeean crash irre y as a 
[1+ Ey — fa " ao 
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where one notes that Bi, = 0 in the expressions for ©, and 6,. The 


] 2 
eigenvalues a, are the roots of equation (17) and initial estimates 
of the roots can be obtained by tabulating the value of the left- 
hand side of Balarion (17) versus a sequence of closely spaced values 
of a5 'S. The values of a, 's at which the sign of the left-hand side 
of equation (17) changes give the upper and lower bounds of the root 
and may serve as a starting point in computing the eigenvalue using 
the variable secant method. 

The limiting case with Bi, = Bi, = Bi (symmetric boundary 
conditions) are of special interest since comparison of the eigen- 
values with the published results can be made. From the symmetry 


condition ¥(O)=¥ (1), One notes from equation (10) that all odd 


terms vanish. Thus, the eigenfunction ¥ <(y) becomes 


(1-0 ./4)*++(4n-3-a./4) 


I) erences aly? )p/ leh (yale) 


After applying the convective boundary condition dy (1) /dy = -BiY;(1), 


one obtains 2, = 051... 


1-a./4)***(4n-3-0./4) a. a. 
| IEE e eget iw en 


“a 
a, we OT ape (2n)! 


The first three eigenvalues are obtained for Bi = 2,20 and the results 
are listed in Table 1 where the resluts from Van Der Does De Bye and 
Schenk [21] and the approximate values of Dennis and Poots [16] are 


also given for comparison. 
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3.2.3 Eigenconstants 
From the orthogonality of the eigenfunctions, the eigen- 


constants can be determined from 


E, = Sgy(1-y)¥ sy )dy/Sqy(1-y¥e (yay (20) 


The integrals appearing in the numerator and denominator can be writ- 


ten as 
, 4ay.(1) 
LoyQ-y)¥5(y)dy = - Pay 
j 
aY. dY. oY. 
ANE eles le ewer oe (2 \a | 
Loy y)Y;(y)dy Za, cre TP ea aye Sa) 4y=0 


Thus, the eigenconstants for the case aA = 0 are 


2 
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The first ten eigenvalues, eigenconstants together with C, AV Ci and 
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The ten eigenvalues are found to give sufficiently accurate solution. 


3.2.4 Location of Point of Incipient Solidification 
Noting that the solidification starts at the location 
where the wall temperature at the upper plate reaches the freezing 


temperature, the liaguid soldification-free length Xe Can be evaluated 
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= EY (1)exp(-a5x¢) (22) 


where 6¢=(T,-T_)/(To-T,,.)- For the determination of Xe it is con- 
venient to define a superheat ratio € as 


e = (Ty-Te)/ (ToT) = 8p. 1 | (23) 


One observes that x, = f(e,Bi 


f 2): 


SS) Analysis in the Freezing Zone 
3.3.1 Formulation and Solution 

The basic assumption that the axial velocity profile re- 
mains parabolic in form even in the presence of axially increasing 
solid layer is well discussed in [1, 7]. In present analysis, only 
steady state solution for Bi,=0 (perfect insulation at bottom plate) 
is sought and no soldification occurs at the bottom plate. The velo- 
city components satisfying the continuity equation dU/dX'+aV/dY=0, 


sou(X' .Y)dY=LU, and the no-slip boundary conditions are found to be 


U(XEY)= =# (Ha -D (24) 
6LU 
VOX Y= BDA - SE 


Substitution equation (24) into energy equation, UdT/dX'+VdT/9dY 
= 0a@T/ay and introducing the dimensionless variables, 
X = X* /(SPPe)* y = Yel; 6 =d/L, o = (T-T.)/(T)-T,) 


one obtains 
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: : . : 2 
a + Yo - 28) - BS (25) 
§ 6 689x496 itay@ oy 


with boundary conditions 
6(x,6)=0 (interface), 36(x,0)/dy=0 (insulation) (26) 
o(O.y)=[(1+e )/e]6 (x¢sy)-(1/e} (initial condition) (27) 


Applying the variable transformation [1, 9] from o(x.y) to 
$(x.n) with n=y/S=Y/S,6 can be eliminated from the boundary condition 
and one obtains the following system of equations which is similar in 


form to the classical Graetz problem. 


2 
page le een 
6 STENT 392 (28) 
6(x,1)=0, 36(x,0)/ay=0 (29) 
(Oon)=[(1+e)/e]o(xesy)-(1/e) (30) 


Using the method of separation of variables by setting 


¢=H(x)N(n), one obtains 
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The solution to equation (31) is 


H(x) = exp [-8°5) 2) 


Hence, the series solution for $(x,n) is 
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al (33) 


where B and N are the eigenvalues and eigenfunctions. 
Using the transformations, 2=8(n-1/2)° and N(n)=w(z)exp(-z/2), 


equation (32) is reduced to Kummer's equation and one obtains . 


wo (1-B:/4)***(4n-3-8./4) 
Ny (m)exoL-8(n-1/2)"/2]-[0 (1+ £ —t yp a(n 1/2)" 


 (3-B./4)**+(4n-1-8./4) 
+ B,/4(n-1/2)(1+ § Cnt By(n-1/2)")] (34) 


J n=] 
From the boundary condition N5(1)=0, D, can be found as 


ie Ae 


/2)%,/%, (35) 
(3-8./4)***(4n-1-8./4) a 


where ©, = 1+ mnti)1 iva : 
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 (1-8./4)***(4n-3-8./4) Ben 
o, = ie ae (7) RU) 


The other boundary condition dN .(0)/dne 0 yields the following equation 


from which the eigenvalues 6. can be evaluated. 
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Thus, one has 


Bt: (1-8./4)+++(4n-3-8./4) 8. . 8B. 


B. (3-8./4)++*(4n-1-8./4) 8. 
ll pe J ajyn 
+ 20,1 ne 2ntT)! a 
RB. 
eer tecn tLe 10 (36) 


The roots representing the eigenvalues can be determined by using the 
procedure described in Section 3.2.2. 


The eigenconstants K, can be evaluated from 


K5 = Lgn(-n)N, (n) @(Osn)an/sgn(1=n)NG (n) dn (37) 


It is noted that o(0,n) depends on the superheat ratio as well as the 
Biot number. The integration in equation (37) is evaluated by a 
Simpson's rule (IBM-SSP-DQSF) employing 200 equal step sizes. The 
first ten eigenvalues, constants D, % Dio and values for dN (1)/dn 


~ dN.,(1)/dn are given in Table 3. Table 4 lists eigenconstants for 


10 
the case Bi, = Q, Bi, =) | peg WON Eee) a1, Ole oe phe 
eigenvalues in freezing zone can be compared with those given by 
= o (constant wall temperature) 


Schenk [15] for the case Bi, = 0, Bi 


] 2 
and the agreement is found to be very good. 


3.3.2 Steady-State Frozen Layer Growth 
Application of the energy balance at the liquid-solid 


interface for steady state yields 


kp 97 (X' 55) / oY = k, oT, (X' 55)/2Y (38) 
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The temperature in the solid-phase layer is described by the following 


system of equations 


Te | dT. (X" 50) 
ioe = OF T A(X me) = 0, Ki eee Vee + ho LT (X sb) =n = 0 (40) 
The solution is 
Bi [(T,-T )/L] 
ate 1+B7 [T-(s/L (y-6) ay 


where Bi. = hoL/k.. 


Using the above result, the steady-state thickness of frozen layer 


5(x) can be found from equation (38) as 


; (14Bi.)ag(x,1)/an 
°C) = 57 Tages 17an-(k,/ Ke) iin 


3.4 Pressure Drop and Heat Transfer Characteristics 
3.4.1 Pressure Drop 

The information on pressure drop is of practical interest 
to prevent the channel from "freeze-shut". The pressure drop is linear 
in the solidification-free zone with a pressure gradient dP/dX 
= (14, )/°.. Defining the dimensionless pressure difference 


(P-P)/ (ou /2), one obtains 
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In order to determine the pressure drop in the freezing zone, 
One integrates the axial momentum equation, USU/3X' + VOU/3Y = -aP/ax' 
+ y(9@u/3ax'* + 3usay’), from Y = 0 to Y = 6 and obtains 


2 ; 
a dy2ay eA cpa yy f52 U dY + yf QUix >6) _ dU(X 0), (44) 
It is now assumed that the axial variation of the solid layer thickness 
is gradual so that Hesyenee (d5/dx')? << 1]. Substituting the expres- 
sion for U into equation (44), one obtains the momentum equation in 


dimensionless form as 


= - 765,04 (45) 
dx os d 39 
Integrating from x = 0 to x and applying the initial condition p 
=144 Pr x at x = O, one obtains 
p(x) = . (T - 1) +144 Pr (x, + 2 = (46) 
6 6 


It is noted that the first term represents the momentum effect and the 


second term denotes the viscous effect. 


3.4.2 Bulk Temperature 
The dimensionless bulk temperature in the solidification- 


free zone is 


Lae 1 
8, = (T.-T)/(Ty-T,,) = Sqledy/J quay (47) 


where U = 6u_y(1-y). Using equations (4) and (5), one has 
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2 
8 = -6 52; E dy. = AO Fit: - 9 
b gar Ey VGC )/dy-exp(-a-x)/ass Oxx<x, (48) 


where dy (1 )/dy = -Bi, YQ). 


For the freezing zone, it can be shown similarly that 


10 a: > ee 
= CX POX aye 
= = - => = >: ] ‘3 i ae ° 
$= (Ty Te)/(Tg-Te) = -6 524 KjdN5(1)/dn-exp[-B5/ 5 : 1/8; (49) 
Noting that > ~ (T,-Te)/(Ty-Te) = (l+e)8, /e-1/e, the expression for 
bulk temperature Sb becomes 
10 ~=dN.(1) ‘=o 
] 6€ Z XOX ,pe 
Ch ee gy eee . = bee 
pb THe ~ Teeghiky an * exPL Big Ps ae ae 30) 


3.4.3 Heat Transfer 
The dimensionless total rate of heat transfer from the 
liquid per unit depth of the channel in the solidification-free zone 


from X = 0 to X is 


Q, = Q/LU PC, (Ty-Te) > OSx<x¢ (51) 


where Q = in - kOT(X,L)/9Y-dX. 


The value of Q, becomes 1 when the liquid is cooled to its freezing 


temperature. Using equation (4), one obtains 


10 aY,(1) 
ee eo ; (52) 


Similarly, the dimensionless heat transfer rate in the freezing zone 
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from x = 0 to x is found as 
Gc 


Thus, (Q, + Q,) represents the total heat transfer rate between the 
thermal entrance (x=0) to any axial location (x=x) in the freezing 


zone. 


3.4.4 Local Nusselt number 
The local heat transfer coefficient, h, is defined by 
Ae kA ak) tht = h(T 


har PO ay O<x<X ¢ (54) 


Spr, 801 (X' 55). : 
q k SY =h(T, Te), X>X ¢ (55) 
Thus, the local Nusselt number, Nu = hL/k, becomes 
Nu = [-90(x,1)/dy]/L8, (x)-0(x,1)]  O<x<x¢ (56) 
Nu =[(1/5)36(x,1)/an]/o,,  x>Xx¢ (57) 


The local Nusselt number provides further insight into the heat trans- 


fer mechanism particularly in the freezing zone. 


e8 Results and Discussion 
3.5.1 Temperature Distributions 

The developing temperature profiles in both the solidifica- 
tion-free and freezing zones for Bi, = 0, Bi, = 2, € = 4 are shown in 
Fig. 2. The difference between the temperatures at the upper and lower 


plates is of interest. The profile at Xe = 0.15 represents the entrance 
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temperature distribution for the freezing zone. For x>0.15, the ice- 


water interface is at the equilibrium freezing temperature T, whereas 


f 


the temperature at the insulated lower plate decreases toward T At 


e 
x = 0.220, the temperature profile becomes fully developed and is re- 
latively flat. 


Numerical results for liquid-solidification-free length x 


f 
versus Superheat ratio € are plotted in Fig. 3 for both the symmetric 
case (Bi, = Bi, = Bi) and the unsymmetric case with Bi, = 0, Bi, = 1, 


2, 10. As expected, Xe decreases with the increase in Biot number 
and increases with the superheat ratio. For a given Biot number, Xe | 
for the unsymmetric case (Bi, = 0) is seen to be about 2 to 3 times 
longer than that of the symmetric case (Bi, = Bi,). 

The axial distributions of the dimensionless bulk temperatures 
0 for Bi, = 0 and € = 1,4 with Bi 


b ] 
Noting that 9 


5 = 1, 2, 10, 100 are shown in Fig. 4. 


rg 8 = 1/(1+€) as xs in equation (50), all the curves 
approach the same asymptotic value depending on €. The development 
length required to reach the freezing temperature o> oe depends on 


Biot number and superheat ratio. 


3.5.2 Steady-State Ice Layer Thickness 

The shape of the steady-state ice layer growth along the 
convectively cooled upper plate is of considerable interest. The axial 
variation of the ice layer thickness for Bi, = 0, Bi, Spalesbec seal With 
fe Fale.) hi 2Ole, 4. (9.15 sand ko /k, = 0.25 is presented in Fig. 5. The ef- 
fects of Biot number and superheat ratio on the ice layer profile can 
be seen clearly. The solution is terminated when § reaches a value of 


3 


10 ~ and the complete "freeze shut" of the channel will not occur be- 
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cause of a constant mass flow rate assumption. With further reduction 
of 6, the assumptions employed in the analysis may not be valid. As 
can be inferred from the characteristic of flat-plate boundary-layer 
heat transfer, the convective heat-transfer coefficient of the warm 
water is high near the solid-phase entrance and the ice layer is seen 
to have zero thickness at x = Xe Further downstream, the convective 
heat transfer decreases with length along the plate and thus the ice 
thickness increases in the downstream direction. The ice-water inter- 
face contour apparently resembles that of forced flow over a flat 


plate [12,13]. 


3.5.3 Pressure Drop and Heat Transfer 

The pressure drop results are given in Fig. 6 for Pr = 13,2 
(water). After the onset of ice formation, the pressure drop increases 
rapidly and the deviation from the curve representing a linear pressure 
drop (p=144 Pr x) occurs. The pressure drop is caused by the flow 
acceleration due to the reduction of the flow area. It is apparent 
that the channel will freeze shut if the pressure gradient is not high 
enough. 

The removal of the sensible heat before the bulk temperature 
approaches the freezing temperature is of practical interest. The 
dimensionless total heat transfer rate Q as a function of x is shown in 
Fig. 7 where the effects of Biot number and superheat ratio are also 
clearly seen. One notes that 0 +1 as 8h > 8 (see Fig. 4), 

Further insight regarding heat transfer mechanism can be 
gained by studying the local Nusselt number behavior. The effects of 


Biot number and superheat ratio on local Nusselt number results are 


aire ie om: 
avis ad 


neg? zt yeyat oot alt ne sonerahe s208q-bh | 
5V rS38vn02 ant “9893 er dade? a = ane “ 
eategds: eur bre soy afd gnols atpnst ixtw 298 | 
~esnt yeten-oot ait ambit st em ato 
tothe rove. woth BStohito tade. satdndat cea 
bac wie “a0 30 
yen tT pr qort wirelen ba ey | 
S.ff «© 19 wots oft eh weetp aie 2d ivze4 gol eruezeny ay — 


eeesoront genb sruecand art <a Weare? oat Yo Jszmo ont NQITAY ae | 
aw2eiy sheath % oa theese ayy silt mort nohistvab sit ‘ne ha if ot 
wol? ott “4 beeuan ah sovb euezary oT vawssa| (s 14 erea) ae 

tus TEN sfoat 5816, Wort sit To wottoubet ont ag. oub ne + inf oi . 
Apid tat 24 tostborp o1waze79 say Fi dure cnet rhtw fanned oly ed 
| aa ~  ewome: 

owseroqnes Atud aid oneited dast gidtease ony 6 Taworiy attri che 
ofT .jearatnt Feohdon te et avutetoqmat pnixsett aie mE : 

nf nwote 2f x te nobronut © 26 5] ater ytenet See, fader: eel 
cele ove off67 conte tai adam Sots to: baits on ™ ' 
(sph soy gels ye en/f =D 20h eaten! an nae a 

ad nba wetnniosm qaameNs teed allele iad nny - | 

¥o eto5tte SAT oF ved “-9gidrruiet Socket <cmienan at 
ows attuans a ane ‘m 


7 y iz ‘ 7 7 7 ‘4 ; 


Shown in Fig. 8. As in the case of classical Graetz problem, the 
local Nusselt number decreases with the axial distance in the soli- 
dification-free zone. The increase of local Nusselt number in the 
freezing zone is related to the increased axial velocity due to the 
presence of the growing ice layer. Apparently, the increase in 
Nusselt number in the freezing zone is related to the corresponding 


increase in the pressure drop. For the cases, Bi, =¢0,eBasts de 2 


2 
and € = 4, the liquid-solidification free zone is long enough so 

that the fully developed Nusselt number is already reached before 

the freezing zone begins. It is noted that the local Nusselt number 

in the solidification-free zone is independent of superheat ratio. 

The local Nusselt number can also be defined using the am- 
bient temperature T, as the reference temperature throughout the whole 
convection cooling section. The numberical results for Nu thus defined 
are plotted in Fig. 9 where the effects of Biot number and superheat 
ratio are seen more distinctively than those shown in Fig. 8 but the 
trend is generally similar. It is of interest to note that the asym- 
ptotic Nusselt number, Nu, = 1.126, in the solidification-free zone 
for! BIL! =2031B1 


] 2 
given by Schenk [15]. 


= 2 agrees well with the asymptotic value Nu, = 1.13 


3.16 Concluding Remarks 

The solution of the steady-state liquid solidification 
problem for laminar forced flow in a parallel-plate channel with 
uniform external convection cooling at the upper plate and perfect 
insulation at the lower plate is obtained by the classical Graetz 


method using the confluent hypergeometric function. The use of the 
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confluent hypergeometric function is particularly effective when the 
boundary conditions are unsymmetric and a large number of eigenvalues 
is required. 

The configuration and the thermal conditions considered in 
the present solidification analysis are encountered in various impor- 
tant engineering applications such as casting of metals in molds, 
heat exchangers using cryogenics, and freezing of ice sheets on north- 
ern rivers, estuaries, waterways and lakes with thermal discharge. It 
must be pointed out that the real physical problems mentioned above 
are of great complexity and the present physical model represents a 
somewhat idealized one. However, the analysis does provide consider- 
able insight into the conditions under which the northern river, for 
example, may freeze completely. In this connection, the pressure-drop 
and heat transfer results presented are useful in estimating the con- 
ditions for the complete blockage of the flow in channels. 


The transient solidification analysis to reach the present 


Steady-state solidification problem is not the subject of present study. 


However, the method described in [8] can be applied to the transient 
solidification analysis involving a parallel-plate channel. 

The developing temperature profiles in Fig. 2 show that the 
flow field is potentially unstable because of top-heavy situtation in 
the case of a horizontal channel. For laminar forced convection in 
horizontal parallel-plate channels, the free convection effects will 
appear only after the onset of a secondary flow in the form of longi- 
tudinal vortex rolls [22,23]. However, the free convection effects in 
the freezing zone are not expected to be important since the tempera- 


ture difference between the lower and upper plates is rather small. 
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The thermal instability analysis is beyond the scope of present in- 
vestigation. 

The experimental investigation for the present problem does 
not appear to have been reported in the past. Thus, the experimental 


investigation is apparently in order. 
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CHAPTER IV 


EXTENDED GRAETZ PROBLEM AND LIQUID SOLIDIFICATION-FREE 
ZONE IN CONVECTIVELY-COOLED PARALLEL-PLATE CHANNELS 


The extended Graetz problem with axial heat conduction is con- 
Sidered for the infinite region consisting of the upstream adiabatic 
region (-© < x < 0) and the downstream convectively-cooled region 
(0 oux o) in parallel-plate channels with fully developed jaminar 
flow. The solution is obtained by the eigenfunction expansion method 
using the confluent hypergeometric function and the Galerkin method 
for eigenconstants. Heat transfer results are obtained for Peclet 
numbers 1, 5, 10, 50 and Biot numbers 1, 2, 10, ~ by using the first 
twenty eigenvalues and the related constants. The axial heat conduc- 
tion effects on the length of liquid solidification-free zone are 
studied and the implications on liquid solidification problem in low 


Peclet number flow regime are pointed out. 
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Greek Letters 


NOMENCLATURE 


parameters in Kummer's equation 

Biot number, hal /k 

eigenconstants, eqs. (3) and (4) 

constants in eq. (10) 

overall heat transfer coefficient between inner wall and 


ambient defined by -k(3T/3Y) h.(T 


fete pO 
thermal conductivity of fluid 


vel 7 ') 
half-distance between parallel plates 

confluent hypergeometric function 

Nusselt number, h(2L)/k 

Peclet number, Ay L/o 

eigenfunctions, eqs. (5) and (6) 

fluid temperature and liquid freezing temperature 
uniform entrance temperature and ambient temperature 
average velocity 

solution of Kummer's equation 

rectangular coordinates 

(X/LPe), (Y/L) 


liquid solidification-free length 


thermal diffusivity 
eigenvalues in eqs. (3) and (4) 
superheat ratio, (T) - PVA ex TS) 


dimensionless temperature difference, (T - T,,)/(Ty - T,,) 


bulk temperature, eqs. (22) and (23) 
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4.1 Introduction 

Theoretical analysis on liquid solidification inside circular 
tubes [1-5] or parallel plates [6] with forced flow is of considerable 
practical interest. The prediction of solid phase formation within the 
tubes or channels is closely related to the well-known Graetz problem 
[1,6]. The analyses [1-6] are all based on the assumption that the 
liquid temperature is uniform at the thermal entrance X = 0. Recent 
theoretical investigations [7-12] on Graetz problem with axial heat 
conduction considering both the adiabatic region (X = -~ ~ 0) and the 


thermal entrance region (X = 0% ~~) show clearly the importance of up- 


Stream heat penetration through X = O for low Peclet number flow regime. 


The analyses [7-12] reported in literature are all limited to the cases 
of either uniform wall temperature or uniform wall heat flux for the 
heated or cooled section. In many practical situations involving a 
cold environment, one may wish to evaluate the effects of the thermal 
insulation and the surrounding sink temperature. Thus, the convective 
boundary condition is important and was considered in [7,13,14] for 
thermal entrance region problem. 

The purpose of this investigation is to study the axial heat 


conduction effects on thermal entrance region heat transfer and liquid 


solidification-free zone in convectively-cooled parallel-plate channels. 


The problem does not appear to have been studied in the past. The in- 
vestigation was motivated by a desire to determine the conditions under 
which the assumption of uniform liquid temperature at the thermal en- 
trance X = 0 is applicable for liquid solidification analysis in paral- 
lel-plate channels. This investigation will shed some light on liquid 


solidification in ducts or channels with forced flow for low Peclet 
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number flow regime where the axial heat conduction is important. 

The elliptic-type energy equation is solved by the method 
of separation of variables and the eigenvalues and eigenfunctions are 
obtained by using the confluent hypergeometric function [15-18]. The 


eigenconstants are determined by using the Galerkin method [19]. 


4.2 Theoretical Analysis 
4.2.1 Problem and Solution 

Consider a fluid with constant physical properties flowing 
between parallel plates from the upstream adiabatic section to the down- 
stream convectively-cooled section (see Fig. 1). It is assumed that the 
flow is a steady plane Poiseuille flow with uniform temperature To at 
X = -2, The Poiseuille velocity is given by U = 3U,[1-(Y/L)°1/2. Ne- 
glecting viscous dissipation, the energy equation with axial heat 


conduction in dimensionless form can be written as 


with the boundary conditions 


68. = 1 at x = -~, 36/dy = 0 at y = O and 1, 


] 


36,/ay =Oaty=0, 36,/ay = - Bie, at y= 1 (2) 


2 
6. = 0. at x =o 


= B55 36,/9x = 36 5/ 9x at x = 0 


where x = X/LPe, y = Y/L, 6= (T - T)/ (5 - T,) and Bi = Not /k 
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and all other symbols are defined in Nomenclature. For the present 
problem, it is shown in Appendix 3 that 6. = 0 (asymptotic solution). 
The solution to the system of eqs. (1) and (2) can be 
obtained by the separation of variables method and the expressions for 
temperature fields, 8, and 855 in the upstream adiabatic and downstream 


entry regions, respectively, are 


2 
; nz BAY by )exp(8a, x/3) -o <x < 0 (3) 


La ge 2 
85 = by CLR (y)exp(-88,x/3) » 0% (4) 


where the eigenfunctions and eigenvalues satisfy the following system 


of equations: 


P+ oOo? - (1 - vA, = 0, ante = 0, nm - (5) 
dy QP M ¥ 
2 
dR 
Deg? jpeoee 2 (a eagle eu (6) 
dy ove, copes eat n 
dR (0) dR (1) 
mod 1 ee SS! eee ee ; 
dy 0, dg BiR (1) (7) 
By introducing the transformations 
z= By’ and w(z) = exp(z/2)R (y) 
equation (6) becomes 
2 B 
DW dw rl _ in G4 c si 
Zz ot (5 < ) dz 4 4 (1 Rj 2 6.) dw 0 (8) 
dz 9Pe 
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which is seen to be Kummer's equation, 


z dw , (b- z) @- aw=0 
dz° dz us (9) 


phereeat= (lave [ice (648¢/9Pe*) (8/4) Ay cay 


The general solution of eq. (9) can be expressed in terms of the con- 
fluent hypergeometric function M(a,b,z) [17]. Thus, the eigenfunctions 


R;(y) are given by 


: eae 2 
R.(y) = exp( BY /2)[C, Mla, s1/2,8.y ) 


1/2 fi 2 
+ Co 58. yM(a, 1/2,3/2,B.y ) (10) 


where C1; and Co are the constants. 


4.2.2 Eigenvalues and Eigenfunctions 

The first twenty eigenvalues and eigenfunctions accurate 
to seven significant figures for eq. (5) with Pe = 1, 5, 10 and 50 are 
obtained in a recent related study [14,20] and will be utilized in this 
study. The first twenty eigenvalues only are listed in Table 1. Equa- 
tion (5) was solved by a fourth-order Runge-Kutta method employing two 
hundred equal steps. 

For the downstream region, the eigenvalues are found by 


application of eq. (7). From eq. (10), one obtains 
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W201 - aay" \Mlay + 1/2,3/2,8,9°) 


2 
“ a.B.yM(a. + 1,3/2.8. 
AC) 5a585Y (a, 1,3/2 BY ) 


+ (4/3) (Cp 589! °y" (a, + 1/2)M(a, + 3/2,5/2,8,9°)] (11) 


From the symmetry condition dk ;(0)/dy = 0, one obtains C 0 (all 


Zita 
odd terms vanish). It is interesting to note that the constants Coy 


correspond to odd eigenfunctions in [12,20] where the upper and lower 
plates are at unequal but constant wall temperatures. Equations (10) 


and (11) now become: 


7 ee 2 
R;(y) = exp( By /2)M(a;s1/2,8,y ) (11) 


ign ee a 2 
dR (y)/dy = exp( Bey ane 8 yM(a;51/2.85y ) 


2 
2. C43.) 3 te? 2 
+ CE 13/2 By )] (12) 


where the constants C1; are being absorbed in the eigenconstants C, 
in eq. (4). 
Substituting eqs. (11) and (12) into the convective boundary condition 


of eq. (7), one obtains 


eee ene \M (an 1/208. M(a. + 1,1/2,8:) = 13 
(Bi Bs 2a, )M(a, iy iz 85) + 2a, (a, ie liye 85) 0 (13) 


where the following recurrence formula is utilized. 
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The first estimates of the roots of eq. (13) can be obtained by tabu- 
lating the value of the left-hand side of eq. (13) versus a sequence 

of closely spaced values of Bs. The values of 8. at which the sign of 
the left-hand side of eq. (13) changes give the upper and lower ends 

of the root and may serve as the starting points for computation. Thus, 
one may inspect the whole spectrum of eigenvalues before the eigenvalues 
are improved by using the variable secant method. | 


The series expansion for M(a,b,z) is 


a(at])***(atn-1) z” 
b(b+1)***(btn-1) n! 


_ a atatl) Zz 
M(a,b,z) = 1+ pee - B(BETY oT . + + e++(]5) 


With, b=, 1/2. a. 


j (1 - AB )/4 and A= 1] + (648°)/(9Pe"), eq. (15) 


becomes 


wo KI-AB.)*** (4n-3-AB . ) 
M(a;.1/2,8.) =" nZ1 ore as(an\ie aoe (16) 


= 


The above expression is used in evaluating Kummer's function M(a,b,z) 
Trance eas 2 where the values of a's are mostly negative in this 
Study. For large ee and Ba Ss the value of Kummer's function eva- 
luated tends to be very large and one cannot maintain the same compu- 
tational accuracy as for smaller ais if the series expansion formula 
from eq. (16) is used. For a. < -2, eq. (16) is used in conjunction 


with the recurrence formula [16]. 
(b - a)M(a - 1, b, z) = aM(a + 1, b, z) - (2a -—.b + z)M(a,b,z) (17) 


The Kummer's function is thus evaluated at two smaller ays utilizing 
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4.2.3 Eigenconstants 
The series expansion coefficients B C. are determined by 


applying the matching conditions at x = 0 given in eq. (2) as 
Ger) Be aby On m REMOTE) 


y= 0 (19) 


With the axial conduction term, the eigenfunctions are nonorthogonal 
and the eigenfunction expansion technique commonly used for the Sturm- 
Liouville system is not applicable. To overcome the difficulty, the 
Galerkin method [19] is used in this study. After multiplying eq. (18) 
by Ms and eq. (19) by ae and integrating from 0 to 1, one obtains the 


following system of equations by taking the first N terms of the series. 


N N 

] ] mr 
nZ] CRY ney oe] BS ova mdY =F ndy (20) 
Oe nacre OF ee te (21) 
n=] n?n Oienen y n=l on'n’ 0 ‘nm y 


where m= 1,...N. The system of simulataneous linear equations for 
N = 20 is solved by using the Gaussian elimination method (IBM-SSP- 
DGELG) with a relative tolerance of ae The numerical integration 
is performed with Simpson's rule. The first twenty eigenvalues, 
eigenconstants B, ~ Boo: Cc, Co along with dR, (1)/dy oY dR 9(1)/dy 
for Bi®="15 °2,°10, ~ and Pe = 1, 5, 10, 50 are listed in Table 2. 
The first matching condition 8, = 6, at x = 0 is checked 


2 
for all combinations of Peclet and Biot numbers, and it is found that 
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three significant figures agree except at points on the upper and lower 
plates where two significant figures agree. The agreement for 36, /9x 


= 96 4/9Xx at x = 0 is generally one order lower than that for 6, = 6 


] ge 
One notes that a discontinuity from zero to a finite value exists for 
(38/2y) 1 


ture, four significant figures agree at x = 0. Thus, one may conclude 


at x = 0 and this causes the discrepancy. For bulk tempera- 


that the Galerkin method is very satisfactory in determining the eigen- 


constants involving nonorthogonal eigenfunctions. 


4.2.4 Bulk Temperature and Local Nusselt Number 
The expressions for bulk temperature and local Nusselt num- 


ber are 


N 
Sale ect | ER Ree ae 2 # 
Ye a! 0 (1-y )8,dy mel aI 9 (1-y") ne] BLY ,exp(80,x/3) dy, axes 0 
(22) 
pom Tel Vey) ee (-BB°x/3)dy, O< x<@ (23) 
(Ai AE U Yo) ney on ne*? n J 31a 
Nu = 2(-28,/8y) 1/8, 
i 2 
= [-2 2, CdR (1)/dy exp(-88,x/3)]/65,, 0< x < (24) 


The integration is performed using a Simpson's rule (IBM-SSP-DQSF) 
using 200 equal step sizes. One notes that Nusselt number in the 


upstream region is zero. 
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4.2-5 Solidification-Free Length 

The prediction of liquid solidification-free length is 
important when a severe cold environment is involved. The solidi- 
fication starts when the wall temperature reaches the freezing tem- 
perature. Beacuse of upstream heat penetration for low Peclet number 
Flow regime, solidification may start in the adiabatic region. Using 
eqs. (3) and (4), one may determine the liquid solidification-free 


length x, from the following equations. 


f 
N 2 
8, (x¢o1) = 2: =] + ne] BAY (1 exp (8a, x-/3) » real aX Bi 0 (25) 
N 2 
85(X¢51) a yates: CAR (1)exp(-88,x_/3) » 0<x,<~ (26) 


where 6, = (T. - T.)/( T, - T,) is the dimensionless freezing tempera- 
ture. For the determination of X es it is convenient to define the 


Superheat ratio € as 


eee (Ta- T.0/( Te a Tole= Ca on | (27) 


It is readily seen that x, = f(e,Bi,Pe). In evaluating X for given 
values of Bi and Pe, the superheat ratio ¢ at which Xe = 0 is deter- 
mined first. For € less than this value, the upstream solution is 


used. Otherwise, the solidification occurs in the downstream region, 


4.3 Results and Discussion 
4.3.1 Developing Temperature Profiles 


The temperature profiles at x = 0 with Biot and Peclet 


the mney | 
Dae sont 
ne =the Toe ont |, 
Saned gator ott 

sada 291989 wot pee si | 
ote At Pa stad ” nid amas ‘en ae 
oer nots eot bt bruett oii oak 
A. ay leds! : 


(2s) 0> x mes ania Savot) 9: 


(as) = gis 0 (0 a Me es = giles 

_— | 

2 - (fa 
‘wraqmes guitestt eeshnotensait ads 2h (oe aT TNS cz - yt) ‘m 


wit onhteb o¢. tagimewnod af df oy to wottenhevevab say 109 
a6 & baka 


atl . eo ; 
cs) | bs Spm. (I STNG oP 
} a | 
navip WT x pot teow tive, ml Ae, ra, a]? * 7* teat noe rr Ai 
tsb at 0 = .x AStAW FB 2 other ssarraque odd .99 find 8 No a 
et noftuloe mesvtequ ‘pil gauley 2iat ogap, 2zet-s 708 chk 

iu vie 

norye nagrd-anwob aod at 210290 oe 20) 


numbers aS parameter are shown in Figs. 2 and 3 where the extent of 
upstream heat penetration can be assessed from the degree of deviation 
from the uniform profile 6(0,y) = 1. It is useful to note that Bi =~ 
is equivalent to a constant temperature boundary condition, and Bi = 0 
represents perfect insulation. The effects of Biot and Peclet numbers 
on the wall temperature 6(0,1) are of particular interest. At Pe = 1, 
the axial heat conduction effect on entrance temperature 6(0,y) is 
considerable. At Pe = 50 and Bi = 1, the profile for 6(0,y) is not 
yet uniform. 

The developing temperature profiles with Pe = 1, 5, 50 in both 
the upstream and downstream regions are shown in Figs. 4 and 5 for Bi 
= 1] and 10, respectively. At Pe = 1, it is seen that the uniform tem- 
perature profile representing the pure axial conduction regime is ap- 
proached at the upstream location x * -0.5 but the uniform entrance 
temperature 9 = 1 is nearly approached only at x = -20. The extent 
of upstream heat penetration is considerable with Pe = 1. At Pe = 50, 
the uniform entrance temperature © = 1 is reached at x = -0.1. The 
developing temperature profiles shown are useful in understanding the 


heat transfer characteristics in the thermal entrance region. 


4.3.2 Heat Transfer Characteristics 

The effect of Peclet number on the axial distribution of 
the bulk temperature is shown in Figs. 6 and 7 for Bi = 1, 10 and Bi 
= 2, », respectively. At Pe = 1, the axial heat conduction effect is 
quite appreciable but at Pe = 50 the effect is rather small. The bulk 
temperature curves for various Peclet numbers cross one another at 


some downstream location before approaching the asymptotic value 85,70. 
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The asymptotic trend is reasonable since the thermal entrance length 
is longer for lower Peclet number. The upstream and downstream deve- 


lopment lengths for approaching the asymptotic values 68, = 1 and 0, 


b 
respectively, depend on both Biot and Peclet numbers. The downstream 
development length increases with the decrease of Biot number but the 
opposite trend is observed for the upstream development length. 

The local Nusselt number results for Bi = 1, 2 and Bi = 10, 
co are shown in Figs. 8 and 9, respectively, with Peclet number as para- 
meter. For a given Peclet number, the Biot number effect on local 
Nusselt number is quite appreciable. Near the thermal entrance x = 0, 
the Nusselt number curve is nearly flat for Pe < 50 excluding the case 
of Bi = ». For the semi-infinite region analysis [14] with axial heat 
conduction effect, Nu > as x > 0 Since the entrance temperature is 
assumed to be uniform. This model is apparently not correct when Pe 
< 50 because of "back" heat diffusion through x = 0. In contrast, Nu 
is finite at x = 0 for the present infinite region analysis. The limit- 
ing case of Bi = ~ (constant wall temperature) is similar to low Peclet 
number mass transfer in laminar flow through circular tubes considered 
TAC 

The asymptotic Nusselt numbers based on the agreement of four 
figures for Nu at the consecutive axial locations are listed in Table 3. 
On this basis, the asymptotic value is already approached at x ~ 5 for 
Pe = 1 and x = 0.2 for Pe = 50. The asymptotic value Nu, = 3.7/3 for 
Bi = ©, Pe = 50 from this study agrees well with Nu. = 3.77 for Bi = », 


Pe = o (no axial conduction) given in [13]. 
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TABLE 3. Asymptotic Nusselt Numbers 


Bi Pe = | 3 10 50 
] ipa tote) Te50 seis) 1.334 
2 2.097 2026 2.000 1.986 

10 3.434 SoU 3a200 Sac 
i 4.003 3.874 Sie oa he) Se 1)s 


4.3.3 Liquid Soldification-Free Length 

The computed results for the solidification-free length Xe 
as a function of the superheat ratio € are shown in Fig. 10 for Bi = 1, 
2, 10 with Peclet number as parameter. At Pe = 1, Bi = 10, for example, 
the solidification occurs entirely in the upstream adiabatic region for 
the range of superheat ratios e = 0 ~ 10 under consideration. On the 
other hand, with Pe = 50, Bi = 1, the solidification occurs mostly in 
the downstream region. With Pe > 10, Bi = 10, one obtains Xe = oe 
Fig. 9 provides some insight into the liquid solidification inside the 


tubes or channels with forced laminar flow in the low Peclet number 


flow regime. 


4.4 Concluding Remarks 

Use of the confluent hypergeometric function for the solu- 
tion of the characteristic eq. (6) involving the infinite region proves 
to be very efficient computationally and requires considerably less com- 
puting time than the Runge-Kutta procedure used in [10, 12]. 

The Galerkin method of obtaining the series expansion coef- 


ficients involving the nonorthagonal eigenfunctions is direct, simple 


98 


ae a: ae 
a ue Te 


- he | 

ut if Faget . porn Ke ‘es 
TT aiete eas a "Cee 
at “<n | | won at 
“bee. f ng “we nee : 
-ers.e a a weer 
ee ete Z *t3.6" 
» if »* ' poate ans 
“Hen J “eavtedoryes' 


» #tphel govt- -ootssattebtlcz ‘ena 107 pais 


ie wo? Of .off ni pwede Gy 2 Orsay 7 a 


wy ars. 


— 


efonaxa 10% .OT = FG Te a9 th 


: : ; = i i » rh ee aceon 
4OT nofeey oitstatbs mepaqu ‘alt wh ¥ 290390 ROFIGOTTT DT 
: : vor) a 0 Bee oe... & 

oft wO .noipessbhene: aetew Of 0 = s’e0tisy Feadieqie 30, 


a} ut ie om etude BH hie PSE of uty .T 4 FS , 02 = at da tW. be 
> 4x ribet We 2 8 oF 5 oT TW nother Maa 
i) shtant setteoreehehle Brat) ort eget tdptent emda saison 
~otitiun tefueS wal ‘atti wh WOT? “sare beon07 Ao hw ataangl 6 


iy rT a 7 


adams gat bil 
fae at} 1% porto’ Seibemdep sql Sraultnos sit Te: ae 
satin not oa ain igt a Benet (a) .ps sheretse 


a o 
or) eae! eben ahiuebs bre! ut fendttstuqms 4 re v 


[St .or} Wr bse srubeoorg ee “ a 
-¥909 notensaxs a an ees Ye! boris a 19° oi 


and more efficient as compared with the rather tedious Gram-Schmidt 
orthonormalization procedure used in [10, 12]. 

Present numerical results show that when the axial heat con- 
duction is included in the energy equation for the extended Graetz 
problem, the assumption of a uniform entrance fluid temperature at 
x = 0 is incompatible with the physical situation and a uniform tem- 
perature must be prescribed at x = --., 

The present analysis does not consider the region x > Xe 
where the prediction of the solid-phase profile is of primary interest 
[1-6]. However, the present investigation does provide some insight 
into the scope of the problem to be expected when solving liquid solidi- 
fication problem with axial heat conduction in a convectively-cooled 
pipe or channel by a numerical method. The heat transfer results pre- 
sented are useful in assessing both the Peclet and Biot numbers effects. 
A discussion on the analysis in the freezing zone (x>x¢) is given in 
Appendix 4 for future reference. 

The Biot number in this analysis corresponds to the wal] 
Nusselt number defined in [13] and the ambient side Nusselt number 
used in [7]. The problem of laminar flow mass transfer in a flat 
channel with permeable walls [18] is analogous to the present heat 
transfer problem with a convective boundary condition. Thus, the 
wall Sherwood number is analogous to the Biot number. The theore- 
tical technique used in this analysis can also be applied to problems 
of heat and mass transfer in a circular pipe. The relevant equations 


for the case of circular pipe are given in Appendix 5. 


ree oe 
+8 wud sieges: ou denser _ unuee a 
“wet movin & bres bah fast eule edt on 1 it ' q m 
= = x te badbiogetg od fit 
4 < = notgey oct Sebtenes aoe beat sede 
Yiolen?. amoz abt iia 23908 aaisaprszovat lade ve 
bh low bitipptt yniiviog aaiw had aegas ad oF" "ae a > a 
bslooo-y slevit 2e4NOD “Bw ful nord ouBMS salen tite = 
ag etlueey ~etanedd phar. ott -hodttan fas tosemun 6d Tannery 
oe zed 11a baw det ae arts dtod gat easezs nt iuts2u'sa ‘ 
nt navtp ay (. xe Si 91105 of sane ont mh eraylens: afi 19. AON; 
Leonena tet abil 10Y te 

rte afd og sbraqeanios shivline atid a “vachimnin, dota why 
‘redinun 3 Taz2utl gbte Inara: ans bis fer}, at bot teb:-v ah om 
tat? ont vena 22st wot? wectast te moTdéorg gfT ieee vale 
‘sot tetszerq of? oF euopplans 2t [Sf] atten oldeemreg Ad hwh rina ‘ | 
| ol? .cudT .aotdbbase yrsenuod svido@vage & at hy eo = | 
-svo%t) oT todd OFA Set 3 epoporians: it: i 
smeTdang: 02 salt 3d oats. neo gout ens 2tdy at bail 
angttaups mevalery out say ‘elutto. & at 1 m 


a eae nevty st sa len, seat ; dhe 


On Oe -_ ‘a a ens oh terh 
i 7 nite io ia, We 


~~? ; 7 r 7 
iP to : mr 
if ur 7 a) 7 : : fail, : 


REFERENCES 


Zerkle, R.D. and Sunderland, J.E., "The effect of liquid solidi- 
fication in a tube upon laminar-flow heat transfer and pressure 


drop", J. Heat Transfer 90C, 1968, pp. 183-190. 


DesRuisseaux, N. and Zerkle, R.D., "Freezing of hydraulic systems", 


Can. J... Chen. Eng. 47, 19695 top. 233-237. 

Zerkle, R.D., "The effect of external thermal insulation on liquid 
solidification in a tube", Proc. 6th Southeastern Seminar on Ther- 
mal Sciences, 1970, pp. 1-19. 

Lock, G.S.H., Freeborn, R.D.J. and Nyren, R.H., "Analysis of ice 
formation in a convectively-cooled pipe", Heat Transfer 1970, Vol. 
1, Elsevier, Amsterdam, Cu 2.9. 

Hwang, G.J. and Sheu, J.P., "Liquid solidification in combined 
hydrodynamic and thermal entrance region of a circular tube", 

Can. J. Chem. Eng. 54, 1976, pp. 66-71. 

Lee, D.G. and Zerkle R.D., "The effect of liquid solidification 

in a parallel plate channel upon laminar-flow heat transfer and 
pressure drop", J. Heat Transfer 91C, 1969, pp. 583-585. 
Schneider, P.J., "Effect of axial fluid conduction on heat trans- 
fer in the entrance regions of parallel plates and tubes", Trans. 
Bivensoc. Mech. sEndcs.-/9.0 toe, pp. '/05-/,/3. 

Agrawal, H.C., "Heat transfer in laminar flow between parallel 
plates of small Peclet numbers", Appl. Sc. Res. A9, 1960, pp. 1/7/- 
189. 

Hennecke, D.K., "Heat transfer by Hagen-Poiseuille flow in the 


thermal development region with axial conduction", Warme - 


100 


-totto bruptt to 
e2e0'79 bas. aes 


pee 


om S , i v 
7 * amet 242 af r aervbust %» > intaaen Pl ont ae paw 
MESES aq ia 03 inet 
. ee eas 7 i 
bruptt po notis Tue iwiniead hsmvsitxs 9 0 Tt oe 
rent 


idl no ashi? o 199 gnortuod asd .2079 aut at notyea t¥ibt fo. 
‘a ot fi 
- OT =F 9g: BRE gona? Ta 
93t to 2tevisoh" , Hoh nell bats Ay eet th 2 a 3 a a 
Tov ,Ov8t Ywetener) 198K y “antag ba fon: ae’, ‘guigoawnes é at pha ay 
: me | 


2 ySrug emerald 0h itveat9 Ww 8 


bankdmoa nt notsagie Fok Pow Ehapt ys te t cusrt2 bis & a ePnewet 
| i 


“edut velac i> & We No! “ps4 Sora y tas ‘emmatt DnG otmeeybonbgt ‘e , 
ie. Ae? laa ; 
ce ee amet <6 | 
A 7 sts 7 


sobyenrtibitee btup it to soatta oT” . 0.4 aTingt brs 2. q) 99d 


oe, 
bas Yetenerd, D600 siatireo niet poqu Tenner Rakion Tat fq ial 
iy! 


B-£62 .qqy ,MIRE.<oTe yotansrt via + "aot 


-ensys teed no nolyouliaty HAT? Tatxe to tagTta" 5.0.4 enne 


Ls 


senest , “eaduy bie espa, Taf | ev8q to anorpsy SORENTO ads at sy 
ANOS 41s vee! BT 203 ‘daa 9 Ae 


el VAG naw? ad wor? ‘vert nh yotensd doa" a8) 7 : Y a 
TT qq 098) 0A , ze 38 fg eda Sot ast fe ge to 298 


4 vee » a 


pact weg, ne at NE 


ots wt wor? of tyabtot nepal yd vatenerd | 
- oil, ‘norsighae> Noh ast wt 


ae 


mo 


il 


Se 


14. 


Leys 


18. 


19; 


10] 


Stoffubertragung, Vol. 1, 1968, pp. 177-184. 

Hsu, C.J., "An exact analysis of low Peclet number thermal entry 
region heat transfer in transversely non-uniform velocity fields", 
AIChE Journal, Vol. 7, 1971, pp. 732-740. 

Deavours, C.A., "An exact solution for the temperature distribu- 
tion in parallel plate Poiseuille flow", J. Heat Transfer 96C, 
1974, pp. 489-495. 

WiUssReo 2; Cheng, K.C. and Ou, J.W., "Low Peclet number heat trans- 
fer in the thermal entrance region of parallel-plate channels with 
unequal wall temperatures", Can. J. Chem. Eng. 54, 1976. 

Van Der Does De Bye, J.A.W. and Schenk, J., " Heat transfer in 
laminary flow between parailel plates", Appl. Sci. Res. 3A, 1952, 
pp. 308-316. 

Hsu, C.J., "Exact solution to entry-region laminar heat transfer 
with axial conduction and the boundary condition of the third 
kind", Chem. Eng. Sci. 23, 1968, pp. 457-468. 

Lauwerier, H.A., "The use of confluent hypergeometric function 

in mathematical physics and the solution of an eigenvalue pro- 
blem", Appl. Sci. Res. 2A, 1950, pp. 184-204. 

PirklesusG. dnd. sigil lato ,eveG... ~A variational approach) to 

low Peclet number heat transfer in laminar flow", J. Comp. Phys. 
9, 1972, pp. 207-221. 

Davis, E.J., “Exact solutions for a class of heat and mass trans- 
fer problems”. Can. J. Chem, Eng> °51, 1973, pp. 562=5/7c. 

Walker, G. and Davies, T., "Mass transfer in laminar flow between 
parallel permeable plates", AIChE Journal 20, 1974, pp. 881-889. 


Kantorovich, L.V. and Krylov, V.I., "Approximate Methods of Higher 


We 


i J my sy 
as We 7 
a fa ei is 


a ee a i K aK 
i 
dns rewibhd vow ists 


."ebfat? Nd tzotov 


i ic iu : ii 
avdiaserb sao ‘orts a? wakgutos oe 
ane 

. Jee ‘eatans"T tool ce wor at itvaet 
“ 
2st dear Tadic tet oed Wes”) MG i thm, ak oa ; 
ddiw 2fannany otal (fe Tags te watpe sone foam 9 a4 


ate . bc and ea A 69 | eveqet thew. fou pam 7 ! 

at “stenend soot “ 5.0 setter bits x “f) 2500 vat naw ef ; 
S30" AS lee . tad (PaGA ,eelate tol ull nemiod wol? venta» _ 
ate-a0¢ AG bo 4 ; 


yo tanwd. taoct thOrmst wol Ba yitne oF nokvolne $o8x3" ‘ b. > ‘Sai $f. 


ij 
oC 


putas afd to notehbaas Yisbntiod oft bas nottoubios fptxe sta tw 

| ROR CIA ae Ade! ES J toe sod nsdd | “baba a 

Collie yr? BmMO seep rienttnoo Yo aaa gat" alt nna | 

-ovq aulevnapls we Te gartehod o67 tne Vee Teahtsredon a ss 
MORSERT .qu ,O88T AS see foe -Tagh me fd 

ad dpeovias lonohSehvay AY, 2.0 ,ottibeate Bas COG iia) 

mil .ghod 40 WOR vente nt satenewy jeer ~adonen tof299 wo 

AS31005 4g NU 

“eney) eso brs teed Be Bent ae 97 anotwytng $56x3" dhe 

SNe BOE .qd) SENOT he. -p0009 mont b. 163 east 

neswiod wot? ventwel Ger vat enet? zea" e+T <20tvsG bas Oy 
OBG-1ES cq -PLeT A forvuot 3491 “eoaala ote 
youpen Yeo ebadtol =f el. y ro ie a i + 

R : 


a 


20. 


(aN 


102 


Analysis", Interscience Publishers, 1958, pp. 54-56. 

Wu, R.S., Ph. D. Thesis, 1976, Dept. Mechanical Engineering, 
University of Alberta, Edmonton, Alberta, Canada. 

Tan, C.W. and Hsu, C.J., “Low Peclet number mass transfer in 
laminar flow through circular tubes", Int. J. Heat Mass Transfer 


5s GAZ sepp. 2187-2201: 


oer 2eeM dull ‘th 


* 


orn aa fic 2, ix 


pie ube 
a 
P > 
i 
\ 
’ a . 


i an vw AUT. a 


7 ws ke - 9 7 ous fl can 
Pe er i, os 


to : 
77 ih - 


ae 


E 
aah " 
= ao = 4 

\! aay 


a tl std re 


if 


=| 


Ww On DW & Ww ld a 


19 
11 
12 
13 
414 
15 
16 


Table 1. 


Pe = ] 


G0. 39604908+06 
0.170589 2E+01 
0. 1550091E+01 
0. 18923958401 
G,2181588E+01 
0. 2436687E+01 
Jo 26674998401 
J. 237985978401 
Je 3977619E+91 
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O.3439175E+01 
0. 369603865+01 
0.3766180E+91 
9. 397194655+01 
O,40669778+91 
G. 42093228401 


17 0.43470098+01 


18 
19 
20 


0.44894678+91 
0.4610063E+01 
0.4736114H401 


Eigenvalues for Adiabatic Region 


5 


9.1598323Et91 
0.2678328E+01 
0.3696697Et91 
OeU345450R401 
0.4976333E+01 
0.55361838+01 
Ve 6U44596E+01 
§9.65133508+01 
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0. fIS 277 Jat 
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0.946 2400E8+01 
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0. 1006569E+02 
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0. 10634877E+02 


10 
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0.53605803 E+01 
0.6355316E+01 
OO, 12407486840) 
0. 79889978401 
O.8692240E+4G1 
0.9344 6008E+01 
3.99535C4 E+GO1 
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9.710728 7E+02 
S. 11592178432 
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6,13026908+02 
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Oo. 39C1528402 
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GW. IES 2ZES5EtO2Z 
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50 


G.7051927E+02 
0.158414 7E+02 
©.1786702E+02 
CO, 1885437E+02 
C.1982894E+02 
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0.23474190E+4+92 
0. 2463882E+92 
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9.2891733E+02 
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0.3085812E+92 
0.33 787055402 
0.3269101E+02 
Ue 2B DUad 8 Ete 
G.3443091E+902 
0.3526995E+02 
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Table 2. 


=], Pe=1 


By 
2.5244216E+00 
0, 15S 4972401 
0. 15393308401 
O.18782668+01 
0.216706 18+01 
0.24223778401 
6. 26536668+01 
5. 2866555E+01 
CG, 30648 1GE+01 
0.3251072E+91 
GO. 34272755+01 
9. 359488658+91 
Os 37IS5087E+0 1 
G.399867CE+01 
0. 4056492E+01 
6.49199125E+01 
Sel SSTCTLETOA 
0.4470782E+01 
©. 4600608E+01 
5247268 74E+01 


Bi=1, Pe=5 


Wcmn DW SW po - 


3, 8787878E+00 
ve 234577 2Et+01 
O,.3318760E+01 
C.49944Z6E+904 
8. 4752435E+01 
©. 5332484E+01 
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0. 6785803E+0 1 
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3. 965COT1ETO1 
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20 6. 10525538402 


Eiqenvalues and Related Constants for 
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dR (1)/dy 
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395 100 29 742701 
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—-9. 7190 1867E+61 
0.190C01760E+01 
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0.163 0799E+01 
-0.1059364Ft01 
0. 1049576E+01 
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Oe TPTI2 TE +6 1 
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0,.1008574E+01 
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0.1007694E+0 1 


1, 2, 10, ~ with Pe = 


B 
n 
-0.7524885£E+00 
9.1554215E-01 
-0.44929 1452-02 
G.2078580E-02 
“0. 12040 935E-e 
0. 786 7793E-93 
= 0333 70GE-U3 
ve 419059 1E-03 
“93708507 EUs 
0. 2519028E-03 
ae 205015 2E=02 
GO. 1701994E=-03 
=9 5 14352675-03 
9.1226706E-03 
~G.1060545E-03 
9.9269096E-04 
-0,8155575E-04 
0.7237496E-04 
-0.6466 1555-04 
0.5812442E-04 


-9.2641423E+00 
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PF sITSSTZOR-O2 
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=e Toe Z E-Ue 
Gat ISIGSSE-02 
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—0..46077574E-03 
VeStal 22 25-G3 
=0..3016407E-03 
0.2634872E-03 
“Oe2321386E-03 
6.2060779E-03 
=). W84Z278C0E-03 
0.1665507E-03 
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0.2776517E+09 
=95 220150 18=01 


C,6991974E-02. 


=G§2 2693922E-02 
0.14964917E-02 
-9.94718716E-03 
0.6518684E-03 
-0.4754652E-03 
0.3618933E-03 
-C.2845650E-03 
0.2295891E-03 
=O 18 9TO07SE-USs 
VU. 1584657E-03 
=0. 1387182E-93 
CeTISoagsE-o3 
-0. 1008691E-93 
0.8858623E-04 
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0.7979801E-04 
-0.6457625E-04 


6.8170123E+00 
=O 3623 71960E-01 
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=F e FIS3S6C3E-— G92 
Cel S13393E-02 
=O 5025585 1-03 
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Oe sS5 153 25-03 
Vane 2 T9259 8-03 
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“Hage C23 rE-03 
9.2040766E-03 
-0.2047998E-03 
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Table 2 continued 
Bi=1, Pe=10 


n Be 


1 0.95885U9E+00 -0.61C 0638E+00 


2 0.3084710E+01 
3 5.4491519E+01 
y 6.5672634E+91 
5) 9.656719052+01 
6 0.73965135+01 
7 &.8148843E+01 
8 ¥.8838795E+01 
9 0.9479480E+01 
10 6.1008003E+02 
971. 9.10647108+92 
92) 0. PHISS68E+02 
43) 0. HPGSSE5E+0 2 
14 G.12192105+02 
15. 9.7266550E+902 
16, 2. B3IZISTEt02 
V7] &. 135630258402 
178) 9. 13990278+02 
419 9. 14404388E+02 
20 9. 1480793E+02 


Bi=1, Pe=50 


41 5.9981068E+90 
2 0.4429122E+01 
37 0. 74856342540 1 
“# 9.1006235E+62 
5 06.1236085E+02 
6 0. 14292658E+02 
7 0.1609381TE+02 
B &.1774397E+02 
9 3.1927251E+02 
10 9. 29069950E+02 
11 0.2294275E+02 
> 02, 23371433 E+O2 
13 0.24524C9E+02 
44 0. 2567992E+02 
15 0. 2678820E+02 
16 3-2 27854178E+02 
47 2. 268822 28+02 
18 9.298/7599E+62 
19 0.3083863E+02 
29 G.3177291E+02 
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dR (1) /dy 


0.1069758F+01 
-0.1103029E+¢1 
0.179340 70E+91 
-0,1067366E+01 
0. 185555325401 
0.10470808+01 
OC. 1040785E+101 
-9.1035946E+01 
0.10321208+01 
—), 10290225490 1 
0.10 2646558+901 
—%.. 102432905461 
0. 170228968+91 
9. 2092328461 
).1019557E+01 
~0.10183565491 
9.1017295E+01 
-90.1710163495+01 
0.107558 12361 


-0.6066959E+G0 
0,10531748+01 
-9.1244283E+01 
0.1303495E+91 
=O. I29 119 Fend) 
0.1258689E+01 
-0.1226433E+01 
0.719916 TEtO 1 
0. 41769138201 
021158778 E+0 1 
0.114 38468+01 
G.323 13988401 
-0.7120832E+01 
0.174119012+07 
-0.1104158E+01 
0.10974 70E+01 
~0.1091278E+01 
GO. 1084962E+G1 
0. 1084902E+91 
0. 10791768401 


B 
n 

—-0.9679815E-01 
0,3317524E-074 
0.7994 728E-02 
0.17657131E-G1 
0.38678 C6E-02 
2.257546 3E-02 
6. 18325558-02 
ie 23 12:18 58-32 
9. 1063004E-02 
0.84872638E-03 
3.69330C0E-03 
0.57698318-G63 
0.4U876694E-93 
0.4176093E-03 
-~0.36162028E-63 
0.31619958-93 
-G0.2788834E-03 
0.2480179E-03 
0.2229359E-03 
0.2052633E-03 


=O, 23555498-302 
O.2492789E-03 
me (2a 83Z23E-02 
©. 26209188-62 
=0.2698127E-G2 
6.1988914E-G62 
-0.1519342E-02 
O, 419696 /E-02 
=0,.966149468-023 
0.7958934E-93 
=0,.6677596E-03 
0.568093 2E-03 
=O AGDIESIE KOS 
0.4317840E-93 
-0.3886232E-93 
0.3605911E-03 
=0.3532824E-03 
0.3711478E-93 
-9,4086 764E-03 
C.4096598E-03 
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Cn 


C.9886881E+00 
-6.7676043E-01 
0. 29598956E-01 
-0.9523924E-02 
G.585069 12-92 
-0.3474181E-02 
©.2382968E-02 
Be. VAZIS I3E=G2 
0.1310616E-02 
=0. 40267 182-02 
C.8257740E-03 
-0,.6784669E-03 
0.5673445E-93 
=0.48151398-03 
0.4139075E-03 
-0.3598028E-03 
0.3160028E-03 
-0.2854661E-03 
6.2502431E-93 
-0.28605005E-03 


0.1083281E+01 
-0.1080273E+00 
0.3520996E-01 
-~0, 1633908E-01 
0.8954873E-02 
=O. 547557 18-62 
0.36339079E-02 
-0.2564836E-02 
9.1898833E-02 
-0.1458775E-02 
O.1154124E-02 
=0. 83509 135-03 
G.7726224E-03 
-0.6489935E-03 
o. 5S22507E-03 
=).47789235-03 
0.4062569E-03 
—0.4159654E-03 
0.1671461E-903 
-0.5221364E-03 
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Table 2 continued 
Bi=2, Pe=] 


n BY 
4 0. 5942040 7E+66 
y) psa oe ahd 
3 Ye 155662258701 
Gg 2. 18883517E+01 
S B. 21737938404 
6 >,2427232E+01 
7 Se 2Z6b5738CE+O1 
8 J. 28695125+01 
9 G6. 3067235E+01 
10 0. 325379068101 
11 ©. 34929013E+01 
12 Ss 80963055704 
493 &. 37568738404 
44 0.390985484+01 
15 0.40575435+01 
16 GO. 420C072E+01 
17 Gs 42379388401 
18 0.447156758+01 
19 6.46013295+01 
20 ©.47275398+01 


Bi=2, Pe=5 


1 0. 1035311E+01 
2 0.2429733E+01 
3 0.3360409F+01 
u 0.41184 70E+01 
5 
6 
7 


eo k 


4768254 E+4+01 
>, 53438275401 
3. 28654055404 
8 9.26344 865E+01 
9 9.6791382E+01 
10 9.7210690E+01 
410. 7697176E+01 
a2 4 279841748 +01 
13 D. 8344284 E+601 
14 0.86895768+01 
$5 we 9O2ZTISSE+09 
16 J+ 93421425401 
17 0.96051964E+91 
46 9.99527T79E+01 
19 0. 1024363E+02 
20 5. 10527032402 


dR (1)/dy 


-0.9265640E+C9 
ePaper ke 4 
0.193 52378401 
eAQISS59EtG 4 
-O.1989438E+G1 
C.19955317E+01 
—-0.1799 853 7E+G1 
0. 200013485461 
“9. 20010368101 
0.200 15578+61 
-9.2001861E4+01 
3.206 2033Et31 
=9, 20024924270 | 
©.20921658+01 
=-9,2002 172E+901 
0.2002158 +61 
-9,20902130Et931 
0.206 20958+61 
~90.200 20538401 
9.2590 2G08E+01 


-0.884U4U404R+C 9 
0.1824 403E+01 
»20138968+01 
0.2041562E+01 
=O. 200320784 7 
C.2Z04C234E+01 
=O sZ2US55S2Ety 1 
0.203 3097E+01 
~0.2030102E+91 
UseZ20Z27S3¢EtU1 
~O e292 5S 25eTU 
0.202 34248+01 
=0. 20277758461 
0.2G020334E+91 
~0.20190665E+01 
0.201798 48+01 
~0.29 1694 2E+01 
0.Z2016045E+01 
~0.,2075239E+01 
0.20714592E+01 


B 
n 


~0.89018682E+00 
2.1983 3718-04 
On DES 182-07 
0. 2814 536E-G2 
v. 765 V659E-02 
ene hee 


war 


by Sp > oot E+ 08 
0.4436 3662-03 
0. 3537323E-03 
On 286679 18-03 
C.2400837E-03 
05 2uze22 oe =OS 
Us A7 B02 268-Us 
6.1593141E-03 


0.13148107E-93 - 


0. 1158693E-03 
279 293768=03 
0.92 C7O3GE-O04 

~8291077E-04 


0.3430053E+00 
0.5661982E-01 


=C374/67 858-04 


a 


2850/7 198E=GzZ 
U, 565926. 7E=Ve2 
Us 23H 1S89E~G2 
e237 268 28-02 
Us a1 Z3482E-C2 
0.73745548E-02 
Us ICO T S458 =92 
0.8964170E-03 
Gs /46T1T3 1E=03 
056307 318R=03 
Js 540223 2E-03 
0.4679 138E-903 
§O92292E-93 
» 3609583E-03 
2320 eee 
-0233733175- 


C 
n 


0,2355142E+06 


—05 27 189S9E—-04 | 


0.8781252E-02 
-0.4013791E-02 
0. Z2249949E-02 
-0,1427062F-02 
9.9819252E-03 
-0.7154U890E-93 
0.5439171E-63 
-G.4271694 8-03 
C. 34423958-93 
-0.2832704E-03 
G.2371704E-03 
-0.2014918E-03 
C.17333538-03 
~150748958-03 
€.1323928E-03 
-0.1173432E-03 
0.1050736E-03 
-0.9696313E-04 


0.7664347E+00 
-C.8461467E-01 
eae ees 
420679 798-01 
‘s TUSS2268-G2 
-C.4554065E-02 
Gs37262625-02 
-O.2273684E-92 
0, 1725836E=02 
=GedSBSISSIE-O2 
0.1989828E-02 
-0.8960718E-03 
0.7497454E-03 
-0.6366286E-03 
0.5474 706E-93 
-C.4769721E-03 
0.4182156E-03 
=0.3711446E-03 
Cx 3336 255E-093 


G.2607667E- aa = 33539078=03 
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Table 2 continued 
Bi=2, Pe=10 


n Bh 


06 TIS2744E+0 1 
Ce 320418 9E+04 
Je 4554698E+0 1 
22 5O4GT46E+O1 
0.6585693E+01 
Je 7T4134176E+014 
Ge 8161599E+01 
Je 834884 3Et+01 
3o 948765 3E+01 
10 %.170086848+902 
11 3. 1065288E+02 
12 5797906 78+02 
13 G.11704802E+02 
14 9.1219596E+02 
1D Pet ZOGoISETU2 
16 04 137250 TEFO 

47 Gs1350583Et02 
18 9. 1399283E+02 
19 5. T4490723E+02 
20 ©. 1481609E+02 


OOon oD UW £& WwW fh = 


Bi=2, Pe=50 


G.1216892E+01 
59.4634 373E+01 
O.7635944E+01 
oO. 1076738E+92 
921237 364E+92 
Ge 1434431F+02 
0.161320Q04E+02 
Se 34777335Et+0 2 
22 1929536Et02 
G6 2071857E+6 2 
11 0.2205866E+02 
412 0.42332785E+02 
93 0. 265357658402 
14-9 ¢25690725+02 
15 Gs 26797208 F02 
16 90. 2786220E+02 
17 0.2888943E+02 
18 05 29882518402 
19 09.3084457E+902 
20 Us 3177833EFO2 


~<A 
QO DYDD UF WD wt 


dR (1)/dy 


-0.8699626E+09 
G.1837647E+91 
<0 2092S FIIETIN 
G22118419E+014 
=“GeZ TOS 27820) 
0.20940973E+61 
“Oe 206337 7201 
O0.2073879E+01 
—-0. 206614 1E+0 1 
0.205978 1E+01 
-O0.2054490E+04 
9.2050031E+O1 
=O. 2046 232E+0 1 
0. Z2042956E+061 
-0.2C40 109F+01 
0.203 7610Et01 
-0.2035400Et01 
3.2033436E+01 
Us 20845 75290) 
0.203 0088E+01 


0.8615370E8+C0 
6.17402719E+4+0 4 
0.2227868E+91 
0.24545158+01 
0.256 0024E+01 
0.24714808+01 
G.2425363E+01 
0.2380583E+01 
0.2341652E+01 
G.2308772F+01 
0.2281099E+01 
0.2257683E+)1 
0. 2237788501 
0.222957 18401 
0.2205532E+01 
0.2192366E+01 
0.218G821E+01 
05277 35768+01 
-90.2160303E+01 
0.2154883E+01 


B 
n 


-0. 1406812E+0C 
0. 5046375E-01 
0.124604 3E-91 

-0.1644055E-01 

“95609 166 15-02 
O.4122499E-02 

See Zoe PS eR Oe 
Oe22 11 8U9E-G2 

“GO prte (239 Bow?2 
VetS?7 7255-62 

=O 092784 1E-G2 
0. 94036075E-03 

=U 779602 158-935 
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Table 2 continued 
Bi=10, Pe=1 


n By 

3269264 18E+00 
2, oes TE +01 
Oe 1633924E+01 
2. 1948412E+01 
G. 22159438401 
9.246006 3E+01 
J. 2683648Et901 
O.2891089E+01 
J. 30852575404 
2. 326844U8ER+01 
9 3442263E+01 


wud COA OO 0) OO UT E Ly OO and 


42 %.3607977E+01 


3 0, 3766605E+H1 
0.39 18967E+01 
ve 4065739E+01 
©.4207487E+61 
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2e4477ITRUE+OT 
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Bi=10, Pe=5 


0. 12557298401 
9,26531168+01 
9.35357725+01 
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6.1925649E+02 
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dR_(1)/dy 


~0.1384287E+01 
0.3936592E+01 
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G.7158693E+01 
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-0.1824189E-02 
02, 15739 782-02 
“0, 137 38278-02 
0.121712547E-02 
=—%., TO831238-02 
0.9833906E-03 
-0.1010893E-02 
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Table 2 continued 


Bi=10, Pe=10 


= 


Bh 


0. 7433665E+01 
Ue 3598)55E+01 
GO. 4896640E+01 
Ge 5839725E+01 
G.67298778+01 
0. I OZ5123E+0 1 
9.825048 1E+91 
G-8921302E+01 
0.9547971E+01 
10 Os ID TST9ISEFOZ 
32 1069684E+4+02 
3011228978 +02 
Ole VIS PIS EPOD 
30 IZ2Z2594E+O2 
0.126958 2E+02 
2213149278402 
Ui 13587868 +02 
0. 14012968+02 
Ge. 14442571E+92 
§.1482713E+02 
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Bi=10, Pe-50 


0.154566 3E+01 
$.517125202+01 
3.38116373=8+01 
9. 1069637E+02 
1.12746118+02 
>.19648458+02 
>.16378288+92 
. 17974168492 
3.19461328+92 
10 9.20857778+02 
41 &.2217 7095402 
12 0. 23429915+32 
43.9.2462474E+02 
14 9. 25768488492 
45.3.2686583E+92 
46 0. 27924565402 
47 0.2894571E+02 
18 0.2993355E+02 
19 0.3089117E+02 
20 9.31821115+92 
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dR (1)/dy 


“9 T29uSS6Ery | 
Oe 3716410E+901 
-0.5802546E8+01 
0.7256938E+01 
-0.8190195Et0 1 
0.8779734E+01 
=0.9157704Er01 
0.949 6496E+91 
x9. 95750458101 
G.9692416E+01 
“9 -9ITTSTISSETO 1 
6.983 7325E+01 
~0 98382778 EF0 1 
9.99917124E+91 
-0.994 3453E+91 
9.99638768+01 
I Comat 
9.999254528+01 
-—0. 10C0202E+02 
0.100 1070E+02 


i 1 20 VISE 
9.32689 26E+01 
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0.6907335Et01 
~9.83157T108F01 
0.92828 42E+01 
-9.9888682E+01 
0.1925607E+92 
-0. 10459 79E+02 
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C 
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0.9346199E+00 


-0.1378446E+00 


0.7792631E-01 
=JersOOSIUSE~OI 
C.2358690EF-01 
-0.1624276E-01 
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-0, 357594 2E-02 
0.3001180E-02 
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Ow2Z20 255 12-02 
-0.1921002E-02 
0. 1694335E-02 
=O. Io 1276 1E-02 
0. 1360423E-02 
<a po2) 7 oe 02 
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“0.252595 198400 
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-0.6149086E-01 
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-0.4228706E-02 
02307 27328-02 
“OAS 10698-02 
0.2665728E-02 
#9527 961908-02 
0.1140741E-92 
“9 S34 19S2ZE-02 
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Table 2 continued 


Bi=e, Pe=] 


n ae 


0.7287309E+00 
9. 13112508+01 
3. 1782419E+91 
0.,.2019032E+01 
Vy ZLIZZIIE HO 1 
D236 285E+01 
9.27587 76E+91 
0). 2964S 15E4+0)1 
Je 315706 1E+O 1 
9.33 38432E+01 
30 351G444E+91 
0. 36744138401 
0.3831371E+01 
H.39821488+01 
2.41274 20E+91 
3242677508 +571 
Ge4¥403610E+01 
G.4535403E+91 
9.46634728+)1 
20 6.4788118E+91 


WOOw ONE Wh 


Bi=o0, Pe=5 


0. 1337162E+01 
G. 2776457E+01 
5, 36868925401 
0.4413191E+91 
9.50360028+01 
0,5589979E+01 
0, 69938838401 
0.6559240E+01 
9.69937358+01 
0, 7402794E+01 
0.77904 258+91 
0 .8159665E401 
3 0,85129178+01 
0.8852091E4+01 
9.9178745E+91 
2.94941715401 
0.97994522+01 
0. 10995518402 
0. 10383135402 
0.10663002+02 


Oo cow HD UE W fy = 


dR (1)/dy 


-9.1542929E+01 
0,.4680637E+01 
-0.7822316E+01 
0.199649008+02 
—9,. 141056528+02 
0.1724730E+92 
-0.2038892E+02 
0.2353654E+92 
-~0.2667215E+92 
0,.2981375E+02 
-0.,3295536E+02 
3. 360 9696E+t92 
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Figure 10. 


Liquid solidification-free length as a function of superheat 
ratio for Bi=1,2,10 with Pe=1,5,10,50 
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CHAPTER V 


CONCLUSIONS 


5.1 Scope of Results 

The present investigation represents a generalization and 
further development of liquid solidification in a parallel-plate chan- 
nel with a constant uniform wall temperature considered by Lee and 
Zerkle [1]. The problems considered in this study are more general 
and the work of Lee and Zerkle becomes a special case with Bi = ~. 
With convective boundary conditions (Newton's law of cooling), one 
may assess the effect of thermal insulation and the results. can be 
used to predict the insulation required to avoid "freeze shut" of a 
channel exposed to a cold environment. Recent theoretical analysis on 
liquid solidification in combined hydrodynamic and thermal entrance 
region of 2 circular tube with constant wall temperature and the related 
experimental results by Hwang and Sheu [2] confirm that the physical 
model used by Zerkle and Sunderland [3] leads to a good approximate ana- 
lytical solution. 

For both symmetric and unsymmetric liquid solidifications in 
parallel-plate channels with uniform convective boundary conditions 
(boundary condition of third kind), theoretical results are obtained 
for liquid solidification-free length, liquid-solid interface profile, 
pressure drop, bulk temperature, heat transfer rate and local Nusselt 
number for a range of values of Biot number and superheat ratio. 

The Graetz problem with axial heat conduction effects in 


parallel-plate channels considering both the upstream and downstream 
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regions has been extended to the case of convective boundary condition. 
It is found that analytical solution in freezing zone cannot be obtained 


for low Peclet number flow regime and numerical solution must be sought. 


5.2 Conclusions and Significance 

The use of confluent hypergeometric function in finding the 
eigenvalues and eigenfunctions proves to be very efficient computation- 
ally and requires considerably less computing time than say the Runge- 
Kutta method. It is particularly effective when the boundary conditions 
at the upper and lower plates are unsymmetric and a large number of 
eigenvalues are required. 

The Galerkin method of obtaining the series expansion coef- 
ficients involving the nonorthogonal eigenfunctions for Graetz problem 
with axial heat conduction effects is found to be simple, direct and 
considerably more efficient than the Gram-Schmidt orthonormalization 
procedure. 

The effect of external thermal insulation on liquid solidifi- 
cation in parallel-plate channels has been investigated analytically 
for steady-state conditions. The Biot numbers,0O and ~, represent per- 
fect insulation and constant wall temperature, respectively. Noting 
that the solution in freezing zone has been terminated at the dimension- 
less liquid-solid interface thickness 6 = ome One may estimate the 
conditions for “freeze shut" using the pressure drop and heat transfer 
results presented. 

The liquid solidification-free length depends on superheat 
ratio and Biot number. On the other hand, the solid layer thickness 


in the freezing zone depends on thermal conductivity ratio kp/k, and 
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Prandtl number additionally. 

The symmetrical solidification problem in parallel-plate 
channels with convective boundary conditions arises in plate coolers 
or liquid-liquid heat exchangers in which the coolant is at a tem- 
perature below the freezing point of the warm liquid. The unsymme- 
trical solidification problem in parallel-plate channels with convec- 
tive boundary conditions such as perfect insulation at lower plate 
and uniform convective cooling at upper plate arises in important 
applications such as freezing of rivers and continuous casting of 
metals. However, it should be pointed out that the physical models 


considered in this study are considerably simplified ones. 


5.3 Some Suggestions for Future Work 

It is desirable to verify the theoretical results reported 
in this thesis by experimental investigation. It is believed that the 
experimental investigation for unsymmetrical solidification with per- 
fect insulation at one plate and convective cooling at other plate may 
serve aS a good starting point in understanding the freezing of river. 
For the analysis on liquid solidification in parallel-plate channels 
with axial heat conduction effects (low Peclet number flow), one must 
use numerical method. 

Future investigations on freezing of liquids in tubes or 
channels may include the case of turbulent flow. In some applications 
it may be desirable to promote freezing such as the use of a portable 
refrigeration device as a temporary shut-off value for a hydraulic line. 
This possibility may have future application in nuclear reactor field 


since the leakage of flowing liquid such as heavy water may be dangerous. 
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APPENDIX 1] 
An Order of Magnitude Analysis of the Governing Equations 


The channel is placed horizontally so that body force effect is 
negligible. For Cartesian co-ordinates assuming no heat sources, no 
pressure work and negligible viscous energy dissipation, the governing 


equations are 


Continuity: ne + 0 (= + a = 0 (1) 
‘Uo pines cll Sai) 8Sl 

Momentum: ae llores OIRO ae (2) 
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pe OV ee) (a) (3) 
ot 3X oY 0 3 ye nye 

c pel [yoo le ge ee Meee aT , aT) (4) 

Nery: Oot aX aY Mice ee 


First, consider the energy equation. Using the following normalized 
variables, x = X/X.s ¥ = VA Mas, Us UU NG V/N Aelita os (Pe a p/0 0 


p= Bi Das > = 6/8. where 9 =T - Le Equation (4) can be written as 


Multiplying the above equation by 1/80 One obtains 
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Ge at /Ye >> 1, then the unsteady term is insignificant. Therefore, 


ey X 7 2 Z 
(aig Gp) a = ral =) aes 
Cc Cac y emp ie): oy 


Assuming the convection terms to be of the same order of magnitude as 


the conduction terms by setting 
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then the energy equation can be simplified to 
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Next, consider the continuity equation rewritten in the following 


form, 
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But if Vo = (=5) VO = -BpV8, one obtains 


or, in dimensionless form, 


Q 


OP Reni, ob Cee ie 
OX Bp ic OX au oy Bp os oy 


where 8 is the coefficient of thermal expansion. 


Using these relationships in the continuity equation, one obtains 
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Multiplying the above equation by X fo U. yields 
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One notes that the ect aric jena exactly the same as that in the energy 

equation. Therefore, if the unsteady term is insignificant in the 

energy equation, it can also be neglected in the continuity equation. 
For water under atmospheric pressure, one has 80 Se USS 


c max 
the right-hand side of equation (6) can be suppressed, since ~ 


V X 
[u ae (hey a is of order one. 
C 


It is also seen that 


The continuity equation now becomes 
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The x and y - momentum equations in normalized form are given by 
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Multiplying the first equation by YU, and the second equation 


2 ; 
by Maeve One obtains 
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Assuming the pressure term to be of the same order of magnitude as the 

viscous term, then one has P_Y e nu Ree nO) OY Par aaUE eA Y - One 
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notices that the coefficients in the momentum equations can be expressed 


in alternate form as, 
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Again, if one assumes that at JY? >> ] and Pr > 1, the unsteady terms 
in the momentum equations can also be neglected. Thus, the x. and y - 


momentum equations now become 
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where the following relationships are used 
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The above equations will now be applied to the solidification-free zone 
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and the freezing zone. 


A. Solidification-Free Zone 
The flow is hydrodynamically fully developed in this region, 


thus, continuity equation becomes 
dv/oy = 0 orv=Q0 


If the flowing liquid under consideration is water (Pr = 7.88 for 
water at 60°F and 13.35 for ice at BOeRY then the axial heat conduc- 
tion effect can be neglected if the Reynolds number is of erder (10) or 


larger. The momentum and energy equations are then given by 


9p/dx = ijn (14) 
ap/dy = 0 (15) 
udd/9x = 3 o/ ay’ (16) 


For this analysis, if one takes Me = L and U. as the average velocity, 


the other reference quantities can be found accordingly. The results 
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Where F is the flow rate per unit depth of the infinite plate, 
vas 2LU. and Rey is a Reynolds number based on Yor However, in the 


text, an equivalent hydraulic diameter is used. 


B. Freezing Zone 
The freezing zone can be subdivided into two parts, the 


liquid phase and the solid phase. 


1. Liquid Phase 
Because of the formation of frozen layer along the channel, 
the reduction of flow area causes the flow in this region to accelerate. 


Energy equation is thus given by 


The rate of solid layer growth is governed by the following energy 


balance at the liquid-solid interface. 
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Where the subscripts % and s denote properties in liquid and solid 

phases, respecitvely, 

6 is the distance measured from centre line to the liquid-solid inter- 
face, and 

L. is the latent heat of fusion. 


Normalization gives 
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where the dimensionless variables are defined by 
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When the solid layer is growing, the coefficient of the middle term of 


equation (17) cannot be suppressed. Therefore, one obtains 
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If 1/Pe << 1, then the axial heat conduction effect is negligible. 


Thus, the energy equation becomes 


2. Solid Phase 


The energy equation is 


134 


sated Oo sin .exgte re j . . om 9d tonnaa (% 


_ 4 2 es, 
ud. beatteb oe 2 es BY 2 


yetwore ef ¥ apust otf 


2emonad nor aun 


at + notaupe ye" 


: : i: ; Ub an 


135 


- oe . We 2 36s 3 “os 
(reptenr it riyolven 20 bye 
Cc Co ox oy 


If one takes 2 = L and Xo = Rey PrL, one notices that Vol Xe = 1/Pe. 
fo 


Thus, if 1/Pe < 1, axial conduction can be neglected in the solid 


region too. The coefficient in the unsteady term is 
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For Ste << 1, the problem is quasi-steady. If the liquid-solid inter- 


face reaches a steady-state, then 


9. r 
Soee 0 and Ate = 0. 


Therefore, the energy equation for the solid phase reduces to a Laplace 


equation as 
or —>=0 without axial heat conduction. (20) 


The interface heat balance equation simplifies to 
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APPENDIX 2 
Derivation of the Expression for K, in Chapter II 


The determination of the eigenconstants K; from Eq. (31) in 
Chapter II resolves to evaluation of the following integrals. 
age ys 
I, Jo (1-n )N,(n)dn 
1, = fp (1-n4)¥,(n)N,(n)dn 
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Using Eq. (26) and integrating, the first integral becomes 
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Next consider the second integral. After some manipulations using 


Eqs.(4) and (26), the integral I, can be shown to be 
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where the following boundary conditions are used. 
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Following the same procedure, I, can be found and is given by 
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APPENDIX 3 


Determination of Fully Developed Temperature 


in Thermal Entry Region 


From physical consideration of the problem, one expects the fully 
developed temperature profile to be independent of x and is of the fol- 


lowing form 
Seater ae (1) 


Eq. (1) must also satisfy the energy equation, 


and the boundary conditions 
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90(x,1)/dy = - Bid(x,1) (4) 


Substituting Eq. (1) into Eq. (2), one obtains 


or Veet yt, (5) 


From Eq. (3), it is found that C, = 0. Thus, substituting ¥ = C, into 


Eq. (4) yields 
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In order to gain a better understanding of the problem, a heat balance 


from X = -e in the adiabatic region to X = X is performed as follows. 
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The term on the left-hand side is the energy input entering the channel 
at X = -~, The first term on the right-hand side is a heat loss through 
convection to the surrounding, while the second and third terms repre- 
sent the axial heat conduction and energy leaving the channel at X, 
respectively. 

Letting 8, = (T., - dle - T.). y = Y/L, x = X/LPe where 


Pe = 400 UL/ky » Eq. (7) can be transformed to 
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With By = 0, one obtains simply 


12-4 fo Sy (XT) dx (10) 


It can be concluded that all the energy entering the channel will 


eventually be transfered out through convection to the surrounding. 
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APPENDIX 4 


Analysis for Solidification with Axial Heat Conduction 


FROZEN 
LAYER 


Fig. Al Configuration and Co-ordinate System for 
Solidification Region with Axial Conduction 


The solidification region is subdivided into two systems; the 
liguid and the solid. The two systems are coupled through the heat 
balance condition at the interface. This coupling together with the 
fact that the solidification boundary cannot be found in advance, makes 
the present solidification problem very difficult to solve. With the 
axial heat conduction effect, the incipient point of solidification may 
occur in the adiabatic region. However, in this analysis, only cases 
where solidification occurs in the convective cooling region (down- 


stream region) are considered. 


1. Liquid System 

As solid layers are being formed along the upper and lower 
walls of the channel, the fluid accelerates in the axial direction due 
to the constriction of flow area. The axial and transverse velocities 


which satisfy the mass flow rate criteria and also the no-slip boun- 
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dary conditions are given by 


LU 

U(x'.Y) = S21 - (D4) (1) 
EUS, 

VOY) = 3 = YT er (2) 


Substituting Eqs. (1) and (2) into the energy equation, 
berm ax" + VT. /aY =~0(8°T./9x'> © 3°T. /ay-) 
x MG Q Q . 


One obtains 
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Introducing the following dimensionless variables, x = X'/LPe, y = Y/L, 


SS vag die by = On - Te)/(T, - Te), the energy equation can now be writ- 


ten as 
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- 1+ ] 
Pcs. 8 (key) (6) 
where e€ is the superheat ratio. 
The boundary conditions (5) and (6) show that the temperature at the 
liquid-solid interface remains constant at its freezing temperature T 


f 
and the initial condition prescribed at X' = 0 is, in fact, the tem- 


perature solution of. the solidification-free zone at x = Xe Thus, one 
has 
CG eaea'a) SeC Rey) exp(-Beex /3) 

Pe n=1 “non ie 


It can be seen that the unknown solid thickness 4 appears in Eq. 
(5). In order to eliminate 6 from the boundary condition, variable 
transformation as described in Chapters II and III may be used. One 


defines 
n= y/d = Y/6 (7) 


The purpose here is to transform 6(x,y) to 8(x,n) so that n changes 
from 0 to 1 when measured from the centre-line of the channel to the 
liquid-solid boundary. 


Eq. (3) now becomes 
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Assuming that the axial solid-thickness variation is gradual, one 


has 
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with the boundary condition 
30(x,0)/an = 1 
o)(x,1) = 0 


>, (0.n) = (Ite )6(x-on)/e cs 


The energy equation then reduces to 
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It should be noted that if the axial heat conduction effect is 


negligible, the energy equation reduces to that given in Chapter II. 


2. Solid System 


For the solid region, it is convenient to shift the co-or- 


dinate origin to the upper plate of the channel as shown in Fig. Al. 


Thus, the steady-state energy equation for pure conduction is simply 


a Laplace Equation. 
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aT 
Kies 80) She ET AX',0) = Ty (14) 


y) 
igiey finite as X' + (15) 
The interface condition which represents the coupling between the 
liquid and solid systems is given by 

, : S 
kylsy-d_ = Kelayel (16) 
Eq. (12) can be written in dimensionless form as 


where the dimensionless variables, o. = (T_ - AU. - Te)» Vitex (Pet 
vie /ieand 6 = o/l sare used. 
In order to eliminate 6 from the boundary condition, a variable trans- 


formation of the following form may be used. 
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ds_y2 
x! dx! 


<<a, 


the energy equation becomes 


2 2 2 
transtoy pine agreceds 1 6 
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Pemeeexia (1-8)" dx! 9x'3& (1-6) med 


with the boundary conditions 


pee Xe ai) oo) 0 (19) 
oan bean 
Be Rs 0) = - LCR", 0) + 2 (20) 
Cap finite as X' +o (21) 
a¢ $ 
and ica =k [= (22) 
y= dy SESS 


From physical consideration, as x' > ™, one expects that the 
temperature distribution in the solid region is developed and is a 
function of y' only. 

The main difficulty in solving this system of equations, Eqs. (8) 
to (11) and Eqs. (18) to (22) is the coupling condition given by Eq. 
(22) since it requires that the temperature solution in the liquid 
phase must also be known. When the axial heat conduction term is 
retained, the energy equation changes from a parabolic-type to an 
elliptic-type. As terms with 6 and dé/dx' are present in the dif- 


ferential equation, one notices that the separation of variable tech- 


a 
" , a 
(ef) 
( "a r) : is . 7 on : 
r ay ; us 
f M4 is A ae r 7 
(93) [~ + (6 3°*) 6] 9 = © (0 eters 
Ls v Te ie? 7 = | 
f; ‘ | ad Y | 
ai o+ '* 86 . | | atiatt = 


i ay 3 


— : 


/ & 7S) ’ i 
ody Jodi stoaqxe eha .@ + 'X 25 .moTtetebleng feotewlq mova. 


4 


6 2t brs bsoo!sveb 2t aorpsn bf foz af? at notjudtgeth oni 6 
«fi yt an 

(8) .2p2 ,2nctraupa to meteye zit patvloe al ytluortrm nisi atT 
oF vd mevio not? Phaes pnt lquon aig 2? (SS) of (81) 4298, bas 
blunt! att of aotfatez stmtersqmad ody dens eortugey ih 

et mad néttsihans toot Tells erid wert nwond sd vo 

hs oo suyd-oT fedievig es int zepant> mote sips wena a 
~tib ddd at jnozang ov 'xb\Ab bas 6 ittiw. amet aA rr ed if, 
“fined wfdstrsy to nattevqse ong yer? seuithon ono ps Saktasy 


147 


nique employed in Chapters II and III is not applicable for the pre- 
sent problem. Thus, with the inclusion of the axial heat conduction 
term, it becomes extremely difficult to seek an analytical solution. 
One notices that the system of equations obtained after the foregoing 
transformation iS no more easier to solve than the original set of 
equations. For the present problem, numerical methods such as finite 


difference method or finite element method are apparently preferred. 
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APPENDIX 5 
Mathematical Formulation for Poiseuille Pipe Flow 


For the corresponding problem in circular pipes, the governing 


energy equation can be written as 


9 oT aT. 1 oT. ap 
2u 01 ‘ (r/r 9) ] a Gl ew vee Nee ik (i mals 2) (1) 


where ro is the radius of the pipe, U, the mean velocity and i = 1] 
Hee euecer LO tnesregi Ons — 25 xe 0) (upstream) and 0 eX 
(downstream), respectively. The appropriate boundary conditions to 
be satisfied are, 
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At n = 0 


6, (0,8) = 6,(0,6) 


Teaser) AS | romney 02D (8) 


Using the same procedure as for the parallel-plate channels, the 


fully developed temperature is found to be 


Spel cee (9) 


The temperature solutions 0, and 85 are now Sought in the following 


forms. 
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On substituting Eq. (10a) into Eqs. (5) and (6), the following charac- 


teristic equation is obtained, 
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Similarly, substituting Eq. (10b) into Eqs. (5) and (7) yields the 
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following characteristic equation 
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Eq. (11) for the upstream region with boundary condition (12) has been 
solved by Hsu [10]. For the downstream region, Kummer's transform can 


again be employed. Thus, one has 
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(£7) 


(eht) 


" rbot alae (Ff 
‘oS miptzagis 2 var ,qorget Meeweannoe gt 103 [ot] ual wl 
eof eno <2udT! seared 


oot ged (ST) mol) i hepogetbetiod Wt tw no | cage 


i 


is SN a Wel 
tO a ae 
. 
SVR ae 4 
ay) (3) 9°. (3),8-90) pea, 
\"9,8- Si | 
>) ets ' 9 + iy P 


7. mG 
5 bo » SVB 
((xYW( gnaw bal slid 9 >" 


enbssdo emo ot 3 olnt zevisavi vel) owt ponte 


pi: 


152 
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which is seen to be Kummer's equation with 


With the symmetric boundary condition dR, (0)/dé = Q, the eigenfunctions, 


then, are given by 
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Differentiating Eq. (17) with respect to € yields 
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The eigenvalues are then obtained as roots of the following equation 


eget ee 
2-O Ole (Pea an One (lt e5/PReq) 


g. [2 2§—n(_* 7" — , 2, 8.) 


2-8 (148-/Pe") : 
Beate (eeepc ste i a (19) 


It should be noted that the above analysis can also be applied to 


the case without axial conduction and with different boundary conditions. 
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DECK FCR COMPUTING EIGENVALUES AND EIGENCONSTANTS 
IN SOLIDIFICATION FREE ZONE 


IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION F(3) 

pO 100 riI=1,5 

READ (5,901) BI 

WRITE (6,1001) BI 

N1=150 


2U EIGENVALUES 


DO 160 JJ=1,20 
READ (5,902) XL 
READ (5,902) XR 


SECANT METHCD 


CALL FUNCT (BI,XL,N1, FL, F,R,N2) 
CALL FUNCT (BI,XR,N1,FR,F,R8, N2) 


30 ITERATIONS 


DO 20 KK=1,30 

XM= (XR¥PL-XL*FR) / (FL-F RB) 

CALL FUNCT (BI,XM,N1,FM,F,R,N2) 
WRITE (6,1004) KK,FM,XM 

IF (DABS (FM) .LT.1.D-7) GO TO 21 
XL=XR 

FL=FR 

XR=KM 

FR=FM 

IF (DABS (FL-FR).LT.1.D-19) GO TO 21 
CONTINUE 

WRITE (6,1005) I1,XM,N2 

WRITE (6,1002) R 

E=XM 


COMPUTE EIGENCONSTANTS 


N=1 

Y=RMULT (E,N) 
REF1=Y 

DO 30 N=2,200 
Y=Y+RMULT (E, ¥) 
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IF (DABS (Y-REF1).LT.0.1D-7) GO TO 31 


39 REF1=¥ 
WRITE (6,444) 

31 Y= (1.D0+Y) /DEXP(E/2.D0) 
IF (E.EQ.1.D0) GO TO 55 
N=1 


Y1I=RMULT (E,N)* (-0.5D0+1. DO/E-SS (E,N)) 

Y2=RMULT (E,N) * (-3.DO0+E/2.D0+2. DO/E-SS (E,N) * (2. DO-£)) 
REF2=Y1 

REF3=Y2 

DO 40 N=2,200 

Y1=Y1+RMULT (E,N) * (-0.5D04+N/E-SS (E,N) ) 

IF (DABS (Y1-REF2).LT.0.1D-7) GO TO 41 


40 REF2=-Y1 
WRITE (6,444) 
41 ¥4=(~0.5D0+¥ 1) /DEXP(E/2. DO) 


DO 50 N=2,200 
Y2=Y2+RMULT (E,N) * (-1. DO+E/2. DO-2. DO¥N+2. DO¥N¥N/E 
X-SS (E,N) * (2. DO*¥N-E) 
IF (DABS (Y2-REF3).LT.0.1D-7) GO TO 51 
50 REF3=Y2 
WRITE (6,444) 
54 Y2=(-1. DO+E/2.D0+Y2) /DEXP (E/2. DO) 
GO TO 71 
cc 
55 Y¥1=0.5D0 
T1=0.5D0 
REF=¥2 
DO 60 NN=2,200 
M=2*NN 
T1=4, DG* (NN- 1) ¥T1/ (M* (M-1)) 
Y1=¥1+T1 
¥ 2=Y2+T1* (M-1) 
IF (DABS (Y2-REF).LT.0.1D-7) GO TO 61 


60 REF=¥ 2 
WRITE (6,444) 
61 Y¥1=(-0.5D9-Y1) /DEXP (E/2.D0) . 
Y2=(-0.5D0-Y2) /DEXP(0.5D9) 
71 EC=-2.D0/ (E¥ (Y1+Y2/BI) ) 


EK=EC*Y 
WRITE (6,222) JJ,EK,EC,Y1,Y2,Y,£ 
WRITE (8,555) E,EC,Y,¥1,Y2 
100 CONTINUE 
222 FORMAT ('0*,*EK*,12,'=", F12.6,10X,*EC=',F12.6,5X, 
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555 
901 
902 
993 
1001 
1002 
1004 
1005 


ce 
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60 
61 


cc 
i 


£*DY=* 7 F 12.60, 5%, *DDY=", F12.6,5X,' Y=", F12.6,5%, "B=", 


XF12. 6) 

FORMAT (! NEED MORE THAN 200 TERMS") 
FORMAT (5D15.8) 

FORMAT (D15. 8) 

FORMAT (D8.4) 

FORMAT (D20. 10,D20.10) 

FORMAT ('1',15X,'BICT NO. =',D10. 3) 
FORMAT (*0',85X,'*R(1) =",D15.8) 


FORMAT (* *,12,5X,'FU =",D15.8,10X,"E =',D15.8) 
FORMAT ('0',80X," EIGENVALUE ',12,'=',F14,8,10X, 


X*n2=* 513) 
STOP 
END 


SUBROUTINE FUNCT (BI,E,N1,FU,F,R,N2) 
IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION F(3) 

A=1,D0 

A1=(1.D0-A¥*E) /4. DO 

NA1=A1 

CALL FCN (BI,A1,E,N1,FN, FU, N2) 
R=FN/DEXP(E/2.D9) 

DO 60 J=1,3 

F (J) =0.D0 

CONTINUE 

RETURN 

END 


SUBROUTINE FCN (BI,T,E,N1,FN,FO,N2) 
IMPLICIT REAL*8 (A-H, 0-2) 

DIMENSION RM1(400) 

FUH=0.D0 

N2=N1 

AM=4,D0*T 

RM 1(1) =AM*E/2.D0 

FN=1.D0+R8M1 (1) 

FCT=2.D0 

DO 10 L=1,N2 

FCT=FCT+2.D0 

AM=AM4+4,D0 

M=L+1 

RM1 (M) =AM*RM1(L) *E/ (FCT* (FCT-1.D9)) 
FN=FN+RM1(M) 
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CONTINUE 
IF ( (DABS (RM1(N2+1))+DABS (RM1(N2))).LT. 1. D-12) 


XGO TO 12 


N2=N2+50 

IF (N2.GE.400) GO TO 11 

GO TO 9 

WRITE (6, 1093) 

STOP 

FU=BI-E 

DO 30 K=1,N2 

FU=FU+RM1(K) *(BI-E+2.D0*K) 
CONTINUE 

FORMAT ("NEED MORE THAN 400 TERMS‘) 
RETURN 

END 


REAL FUNCTION BMULT¥*8 (Z,N) 

IMPLICIT REAL*8 (A-H,O-Z) 

RMULT=1. DO 

DO 37 Is1,8 

M=2*L 

RMULT= (4. DO* 1-3. DO-E) *E* RMULT/ (M* (M-1)) 
RETURN 

END 


REAL FUNCTION SS*8 (E,N) 
IMPLICIT REAL*8 (A-H, 0-2) 
SS=0, DO 

DO 80 J=1,N 

$S=SS+1 ~-DO/(4.D0*JI-3 »DO-B) 
RETURN 

END 
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80 


81 


70 
71 
100 
1174 
222 


555 


CALCULATE "DE*, AND *DNN! 


IMPLICIT REAL*8 (A-H,0-Z) 

DO 100 J=1,20 

READ (5,111) E 

N=1 

C=RMULT (E, N) 

DE=C* (-0.5D0+1.D0/E-SS (E,N) ) 

DNN=C* (2.D0-2) 

REF1=DE 

REF2=DN 

DC 80 N=2,150 
DE=DE+RMOULT (E,N) * (-0.5D04N/E-SS (E,N)) 
IF (DAES (DE~REF1).LT.0.1D-6) GO TO 81 
REF1=DE 

WRITE (6,555) 
DE=(-0.5D0+4+DE) /DEXP (E/2.D0) 

DO 70 N=2,150 
DNN=DNN+RMULT (E,N) * (2. DOXN-E) 

IF (DAES (DNN-REF2).LT.0.1D-6) GO TO 71 
REF2=DNN 

WRITE (6,555) 
DNN=(-E+DNN) /DEXP (E/2. D0) 

WRITE (6,222) J,E,DE,DNN 

CONTINUE 

FORMAT (F10.6) 


FORMAT (*0',5X,*J=",12,10X,*E=",F10.6,10X,'DE=', 
XPF10.6, 10X,*DNN=',F10.6) 


FORMAT (* NEED MORE THAN 150 TERMS ') 
STOP 
END 
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DECK FOR COMPUTING EIGENCONSTANTS, ICE PROFILE, 
PRESSURE DROP, HEAT TRANSFER RATE, BULK 
TEMP., LOCAL NUSSELT NO. IN FREEZING ZONE 


IMPLICIT REAL*8 (A-H,O-Z) 

DIMENSION ESQ(20) ,SQ(20) ,EK1(20) , E2 (20) , EK (20) 
DIMENSION DELTA(15) ,DE(20) , DNN (20) 

DO 3 I=1,15 

READ (5,333) DELTA(I) 

DO 49 I=1,20 

READ (5,444) DE(I) ,DNN (ZI) 


EIGENVALUES FOR FREEZING ZONE CAN BE DETERMINED 
AS IN SOLIDIFICATION FREE ZONE 


DO 52 I=1,20 

READ (5,333) E2(1) 

CONTINUE 

READ (5,111) BI,XX,TS 

READ (5,904) QXF 

WRITE (6,222) BI,XX,TS 

DO 50 I=1,20 

READ (8,902) E1,EC1,Y¥,DY,DDY 

EK1 (I) =EC1¥*Y 

ESQ (I) =E1*E1 

FORMAT (5D15.8) 

DO 100 J=1,20 

SQ (J) =E2 (J) *E2 (J) 

suM=0. DO 

DO 55 K=1,20 
V=8.DO*ESQ (K) *XX/3.D0 
S1=EK1(K) / (DEXP (V) * (ESQ (K) -SQ(J) )) 
IF (DABS(S1).LT.0.1D-6) GO TO 56 
SUM=SUM+S1 

CONTINUE 

WRITE (6,666) K,SUM 

SS= (1. D04TS) *B2 (J) *SUM 

EC2= (2.D0/ (DE (J) *TS) ) *(1.D0/E2 (J) +SS) 
EK (J) =EC2*DNN (J) 

WRITE (6,555) J,EK(J) ,EC2,DE (J) ,DNN (J) 
FORMAT (D15.8) 

CONTINUE 

FORMAT (F9.6,F9.6,F9.6) 
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222 FORMAT ('0',*BIOT NO. =',F9.6,5X,'X=",F8.6,5X, 
K'TS=', F8.6) 

333 FORMAT (F10.6) 

444 FORMAT (F9.6,F12.6) 

555 FORMAT ('0',10X,'J=',12,10X, "EIGENCONSTANT=',F10.6, 
X10X, *EC=",£15.6,5X, "DE=',F 10.6, 5X, 'DNN=',F10.6) 

666 FORMAT (*0',2X,*K=',12,5X,'SUM=',F10.5) 


ce 

ce 

cc K FOR ICE = 1.28BTU/(HR-FT-F) 

a K FOR WATER = 0.332 BTU/(HR-FT-F) 
cc KL/KS APPROX. 0.25 

cc 

cc 

21 READ (5,901) XMAX,PR 


777 FORMAT (F8.6,F10.6) 
IF (XMAX.EQ.-0.D0) GO TO 300 
BIS=BI/4.D0 
TE=4.D0/T 
WRITE (6,888) TE,PR 
888 FORMAT (*'1',10X,'*TW=',P10.6,10X,'PR. NO. =',F10.6) 
I=0 
J=1 
K=1 
IP=0 
IStT=-1 
DEL=1. D0 
DEL1=1.D0 
DEL2=1.D0 
A=0.D0 
B=0. D0 
X=0. D0 
X 1=-DELTA (1) 
PFREE=2.4D1*PR¥XX 


PV¥1=0.0D0 
QOT=0.D0 
4 I=I+1 
TRG(TS1) 5,5,6 
6 J=J+1 
LF S{DELTA (3)) 21,21,51 
51 I=1 
5 IP OCIE) 7,7,8 
7 D=DELTA (J) *{(1.DO0/DEL1) +(1.D0/DEL) ) /2.D0 
IP=1 
GO TO 9 


4 D=DELTA (J) *((1.D0/DEL2) + (4. DO/DEL1) # (1. DO/DEL) ) 
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K/3.3B9 
TpP=0 
3 BN=A+*D 
SUMI=0.D0 
SUMJ=0.D0 
DO 10 M=1,20 
ARG=-8 .DO*SQ( M)*BN/3.DO 
IF (8.DI+ARG) 25,25,22 
22 C=EK(M )*DEXKP{ ARG) 
SUMI=SUMIt+C 
SUMJ=SUNJ+(1C/SQO(M) 3} 


10 CONTINUE 
GO TO 23 
25 WRITE( 6,26) 
26 FORMAT (5X%,'ARG. EXCEEDS 80.08 ) 
23 DN=SUM1I*{(1.D0O+S5815S )/( BIS*( SUMI-TE)) 


TF (DN.LTwj1.D-4) GO TO 21 
ey ( DN-DELt0.159-3) 11,12,13 


i3 IF ( DN-DEL-0.1D-3) 12,12,11 
11 DEL=DN 
K=K+1 
IF (8-15) 14,14,15 
15 WRITE( 6,16) 
16 FORMAT (1X,'NOT CONVERGING? ) 
GO To 21 
14 IF (CIP) 8,8,7 
12 WRITE (6,1000) SUMI,SUMJ 


1000 FORMAT (* ',30X,51426,5X,E14-6) 
QINT=SUMI/DN 
IF (ISTT) 200,200,201 
200 QINT1I=CINT 
ISTT=1 
201 QT=OT—-2.DO*DELTA( J )*( CINTI+FQINT ) 
DELIT=DELI*DEL1I*DEL1 
DNT=DN*DN*DN 
IF (IP) 28,28,29 
29 S=DELTA( J )*((1.D0/DELIT )+( 1.D0/DNT))/2.D0 
GO TO 31 
28 S=DELTA( J )*((1.D0/DEL2T )+{ 4.D0/ DELIT )#(1.2D0/DNT)) 
x/3.2R0 
31 PY=2.-4D1*PR*( BtS ) 
x2=x1 
<i=x 
K=X+DELTA( J) 
K=1 
DPDEL=( DN-DEL2 )/( K-X2 ) 
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20 


13 
300 


D2ZDEL=(( DN-—DEL1 )/{( X—X1 )-( DELI-DEL2 )/( K1-X2) 0/( X-X2) 


PM=6.2D0*{(1.D0—-DELI*DEL1 )/(5.DO* DELI*DEL1 ) 
P=PMt+PV1 

IF (iP) 34,34, 27 

WRITE (6,18) X1,BRBEL1,DDEL,D2DEL,5N 
FORMAT ('0!,2%,F926,4(2XK,E14.6)) 
TB=TS*(1.D0/TS-32.D0*SUNI/2.2.00)/(1..DdD0t+TS) 
TT=1.D0/1S-TB*(1.D0tTS)/TS 
RN=2.D0*SUMI/( DELI*TT) 

WRITE (6,32 )P,°M,PV1,0T1,T8,RN, X 
PORMAT (7( 2X%,F12.5) 

AX=X+XX 

PP=P+PFEREE 

QQQ=QT1It+Cxr 

WRITE (9,905) XK,AX,DEL1,PP,00GQ,TSB,RN 
FORMAT (# 14,8100e3,610.3,5F10.3) 
DEL2=DELI1 

DELI=DEL 

DEL=DN 

DEL2ZT=DELIT 

PY1=PV 

QINTI=QINT 

QT1=QOT 

IF (iP) 20 ,20,193 

A=BN 

B=BtS 

IF ({X¥—-KXMAX) 4,4,21 

CONTINUE 

STOP 

END 
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LIQUID SOLIDIFICATION 
WITH EOTTOCM PLATE INSULATED BI? = 0 


PROGRAMS FOR SOLIDIFICATICN FREE ZCNE 


1. EIGENVALUES 


IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION X(20),F(20) 
BI1=0.50 


DO 200 I=1,4 
READ (5,901) BI2 
WRITE (6,902) BI1,BI2 


EIGENVALUES 
DO 100 J=1,10 


READ (5,901) X(1) 
READ (5,901) X(2) 


SECANT METHOD 


CALL OFCH® (851,812, % (1) ,F (1) , PHL1, PHI2,R1, 82) 
DC 5 KK=1,20 

TI=KK+1 

JJ=KK+2 

CALL FCN (BI1,BI2,xX t{IZ),F (i1)),PHI1, PHI2, 81,82) 
X (JJ) =(X (11) *F (KK) -X (KK) *F (11) ) / (F (KK) ~F (11)) 
DD=X (JJ) -X (IT) 

WRITE (6,904) KK,DD,X (JJ) 

IF (DABS (DD) .LT.0.1D-7) GO TO 6 

CONTINUE 

CALL FCN (B11,BI2,X (JJ) ,F0,PHI1,PHI2,R1, 82) 
C1= (-B11+BI2) * (1.D0+R2) ¥DSOQORT(X (JJ) ) / (PHI1*2. D0) 
WRITE (6,905) J,X(dd),FU 

WRITE (6,906) C1 

CONTINUE 

CONTINUE 

FORMAT (F10.4) 


163 


ne 


sosagre a8 a 
ve =20 — 


S08 atas NOUPADESD 


a : or +e vee 
Ye bated 00S oa 
ett eae 
“TF e20gaveror: | 0 3 
ons"tgbgeeornce 
gi oe 


 tonrE Lb) 


(Sh, FO, SINE FINI, OF) Ty (Fx STENETAD 124 4082 
OS,f=aa 2 Of ~ 

t+xaeIT 

| Seana 

(£6 . Task IRF PERD, (27) %, (Tp RASS eTAy Le) @ 
(tzT} 1+ (aay id dl (At) X= (Aa) a pabenines 


x 08.8) 
Le er a ( 


(om, TR STaa, pres, <0 (66) X, 878.1 TE) 100. 
(00 .O4F THRYA( Ub) 2) PA2 *{SHeOd sf) * (SIash Fe 
UG, alae ry: arr 
a7 TlH 


ee 
ij 


@.0ra 


902 


904 


905 


906 


CC 


10 


20 


30 


903 


cc 
cc 
cc 
cc 


FORMAT ('1',"BOTTOM PLATE BIOT NO. =',F10.6,5X,'UPPER 
XPLATE BICT NO. =',F10.6) 


FORMAT ('0%,"NO. =',12,5X,'DIFFERENCE =',F10.6,10X, 


X*X =',F10.6) 


FORMAT (*0',60X,"EIGENVALUE ( ',12,' ) =',F10.6,10X, 


X'FU =",F12.7) 


FORMAT ('0',60X,'C =",F10.6) 
STOP 
END 


SUBROUTINE FCN (BI1,B12,EE,FU, PHI1,PHI2, 81, R2) 
IMPLICIT REAL*8 (A-H,O-Z) 

DIMENSION 8M1(200),RM2 (200) 

E=EE/4.D0 

RM1(1)=(1.D0-2) *£/2.D0 

RM2(1)=(3.D0-2) *£/6.D0 

R1=BM1(1) 

R2=RM2 (1) 

DO 10 N=2,150 

FCT=2. DO¥*N 

M=N-1 

RM1{(N) = (4. D0*N-3.D0-E) *RM1 (M) XE/ (FCT* (FCT-1.0D0) ) 
RM2 (N) = (4. D0*N-1.D0-E) *RM2 (M) ¥E/ (FCT* (FCT+1.D0) ) 
R1=R1+RM1(N) 

R2=R2+BM2 (N) 

CONTINUE 

IF (DABS (RM1(150)-RM1(149)).GT.0.1D-7) WRITE (6,903) 
$1=0.D0 

$2=0.D0 

DO 20 11=1,150 

$1=S14+RM1 (11) * (EE-8.DO*I1- (BI1+BI2) ) 

CONTINUE 

DO 30 I2=1,150 

S$ 2=S2+RM2 (12) * ( (BI14+BI2) /4.D0-E+2.D0*1I2+1.D0) 
CONTINUE 

PHI1=EE- (BI1+BI2) +391 

PHI2= (BI14+B81I2) /4.D0-E+1.D0+S2 

B= (EI11-B12) *(BI1-BI2)/4.D0 
FU=PHI1*PHI2+B* (1.D0+R1) *(1.D0+R2) 

FORMAT ("NEED MORE THAN 150 TERMS") 

RETURN 

END 


2. EIGENCONSTANTS AND ICE-FREE LENGEH 
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cc 
IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION EG(10) ,ESQ(10) ,¥1(10) ,F1(201) , F2 (201) 
XZ1(201) ,22 (201) 
BI1=0.D0 
DO 200 I=1,4a 
READ (5,901) BI2 
WRITE (6,902) BI1,BI2 
DO 100 J=1,10 
READ (5,901) EE 
ESQ (J) =FEXEE 
E=EE/4.D0 
READ (5,901) C 
NDIM=201 
DO 90 MM=1,NDIM 
K=MM-1 
Y¥=0.05D-1*K 
CALL FN (Y,E,EE,C,W,W1,W2) 
F1(MM)=W1 
F2(™M)=w2 
90 CONTINUE 
Y1(3)=8 
SP=0.05D-1 
CALL DCSF(SP,F1,21,NDIM) 
CALL DOSF (SP, F2,22,NDIM) 
EG (J) =Z1(201) /22 (201) 
WRITE (6,904) J,EE,EG(J) ,Y1 (J) 
100 CONTINUE 
ec 
CC ICE FREE LENGTH 
cc 
WRITE (6,905) 
X=0.D0 
DO 60 1I1=1,50 
X=X+0.1D-1 
T=0.DO 
pO 50 12=1,10 
DOM=ESQ (12) ¥*X 
IF (DOM.GT.8.D1) GO TO 51 
T=T+EG (12) ¥Y1(12) /DEXP (DOM) 
50 CONTINUE 
51 TS=(1.D0-1) /T 
WRITE (6,906) X,TS 
60 CONTINUE 
200 CONTINUE 
901 FORMAT (F10.6) 
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902 FORMAT (*1',*BOTTOM PLATE BIOT NO. =',F10.6,5X, "UPPER 
XPLATE EICT NO. =',F10.6) 
904 FORMAT ("0','J3 =',12,10K, EIGENVALUE =',F10.6,10X, 
X"RIGENCONSTANT =',F10.6,10X,'Y1 =',F12.6) 
905 FORMAT ('1%,! ICE FREE LENGTH’) 
906 FORMAT ('0',5X,'X =',F10.6,10X,'TEMP. OF SUPERHEAT =! 
X,F10.6) 
STOP 
END 
ec 
SUBRCUTINE FN (¥,E,EE,C,W,W1,W2) 
IMPLICIT REAL*8 (A-H,0-Z) 
DIMENSION BN1(200) ,RM2 (200) 
YY=Y-1.D0/2.D0 
RN1(1)=(1.D0-E£) *EE*XYY*YY/2.D0 
RM2 (1) =(3.D0-2) *EE*YY*YY/6. DO 
R1=BM1 (1) 
R2=RM2 (1) 
DC 10 N=2,150 
FCT=2.D0*N 
M=N-1 
RM1(N) =(4.D0*N-3.D0-E) *RM1 (M) EEX YY*YY/ (FCT* (FCT-1.D0) ) 
RM2 (N) = (4.D0*N-1.DO-E) *RM2 (M) XEEXYY*YY/ (FCT* (FCT4+1. D0) ) 
R1=R14+RM1(N) 
R2=R2+BM2 (N) 
10 CONTINUE 
IF (DABS (RM1(150)-RM1(149)).LT.0.1D-7) GO TO 11 
IF (DABS (RM2(150)-RM2(149)).LT.0.1D-7) GO TO 11 
WRITE (6,903) 
STOP 
11 W= (C* (1. D04#R1) #DSORT (EE) *Y¥* (1. D0+#R2) ) 
X/DEXP (FE*YY*YY/2.D0) 
W1=Y¥* (1. D0-Y) *# 


W2=W1* 
903 FORMAT (‘NEED MORE THAN 150 TERMS") 
RETURN 
ERD 
cc 
cc 
cc PROGRAMS FOR FREEZING ZONE 
cc 
cc 1. EIGENVALUES 
cc 
cc 


IMPLICIT REAL*8 (A-H,0-2Z) 
DIMENSICN X(20),F(20) 
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WRITE (6,901) 
DO 100 J=1,10 
READ (5,902) X{1),X(2) 


CC SECANT METHCD 


cc 


cc 


CALL FF (X(1),F(1),81,82) 

DO 5 KK=1,20 

II=KK+1 

JJ=KK+2 

CALL FF (X(II),F(II),&1,R2) 

K (JJ) = (X (11) *F (KK) -X (KK) *F (11) ) /(F (KK) -F (1I)) 
DD=X (JJ) -X (IT) 

WRITE (6,904) KK,DD,X (JJ) 

IF (DAES (DD).L1.0.1D-7) GO TO 6 

CONTINUE 

CALL FF (X(JJ) ,FU,81,R2) 

WRITE (6,905) J,X(JJ),FU 

C=-DSQRT(X (JJ) )* (1. DO+R2) / (2.D0* (1. D04+81) ) 
WRITE (6,906) C 

CONTINUE 

FORMAT ('1°,'FREEZING ZONE*) 

FORMAT (F10.6,F10.6) 


FORMAT ('0*,*NO. =',12,5X, DIFFERENCE =',F10.6,10X, 


FORMAT ('0',60X,*EIGENVALUE (',12,') =",F10.6,16X, 


X'FO =',F10.6) 


FORMAT ('01,60X,'C =',F10.6) 
STOP 
END 


SUBROUTINE FF (EE,FU,R1,R2) 

IMPLICIT REAL*8 (A-H,0-2) 

DIMENSION RM1(200) ,RM2 (200) 

E=EE/4.D0 

RMN1(1)=(1.D0-2£) *£/2.D0 

RM2 (1) =(3.DO-E) *E/6.D0 

R1=RM1 (1) 

R2=RM2 (1) 

DO 10 N=2,150 

FCT=2. DO¥N 

M=N-1 
RMN1(N)=(4.D0*N-3.D0-E) *RM1 (M) *E/ (FCT* (PCT-1.D0) ) 
RM2 (N) = (4.D0O*N-1.DO-E) *RM2 (M) *E/ (FCT* (FCT+1.D0)) 
R1=R1+RM1 (N) 

R2=R24+RM2 (N) 
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CONTINUE 

IF (DABS (RM1(150)-RM1(149)).LT.0.1D-7) GO TO 11 
IF (DAES (RM2 (150) -RM2(149)).LT.0.1D-7) GO TC 11 
WRITE (6,903) 

STOP 

S$1=2.D0*E 

$2=E-1.D0 

DO 20 J=1,150 

S1=S14RM1 (J) * (2. DO*E-4. DO*J) 
$2=S2-RN2 (J) * (-E+2. DO*J+1.D0) 

CONTINUE 

FU=S1* (1.D04+R2)+2.D0* (1.D04+81) *S2 

FORMAT ("NEED MORE THAN 150 TERNS") 

RETURN 

END 


2. EIGENCONSTANTS, ICE THICKNESS PROFILE, 
ERESSORE DROP, EULK TEMP., HEAT 
TRANSFER RATE, LOCAL NUSSELT NO. 


IMPLICIT REAL*8 (A-H,O-Z) 
DIMENSION £1(10),EK1(10) ,£2(10) ,EK2(10) ,C1(10), 


XC2(10) ,DNN (10) ,F1(101) ,F2(101) ,21(101),22(101), 
XEK (10) ,SC(10) , DELTA (15) 


READ (5,901) BI1,BI2 
WRITE (6,904) BI2,BI1 
pO 50 I=1,10 

READ (5,902) E1(I),#K1(ZI),C1(Z) 
READ (5,901) E2(I),C2(Z) 
$Q (I) =E2 (I) *E2 (I) 
CONTINUE 

po S2 I=, 10 

READ (5,222) DNN(Z) 

pO 3 I=1,15 

READ (5,222) DELTA (I) 

DO 300 NN=1,3 

READ (5,901) T,XE 

READ (5,222) QXF 

WRITE (6,906) T,XE 

DO 900 J=1,10 


NDIM=101 
DO 100 MM=1,NDIiIM 
K=MM-1 


Y=0.10D-1*K 
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S=0.D0 

DO 90 L=1,10 

CALL FEF (Y,E1(L),C1(L) ,W1) 
S1=E1(L) *E1(L) *XE 

IF (S1.GT.8.D1) GO TO 91 
S=S+EK1(L) *W¥1/DEXP (S1) 
CONTINUE 

SUM=(1.D04T) *S/T-1.D0/T7 
CALL FF (Y,22 (0) ,C2 (3) ,W2) 
F=Y* (1.D0-Y) *W2 
F1(%M)=F*soM 

F2 (MM) =F¥*H2 

CONTINUE 


CC NUMERICAL INTEGRATION 


cc 


900 
222 
901 
902 
903 


904 


905 
906 


cc 
ft 
cc 
cc 
ce 
cc 
cc 
21 


SP=0.1D-1 

CALL DQOSF (SP,F1,21,NDIM) 

CALL DOSF (SP,F2,22,NDIM) 

EK2 (J) =Z1(101) /22 (101) 

EK (J) =DNN (J) *EK2 (J) 

WRITE (6,903) J,E2(J),EK2(J) ,DNN(J) ,EK (J) 
WRITE (6,905) 21(101),Z2(101) 

CONTINUE 

FORMAT (F10.6) 

FORMAT (F10.6,F10.6) 

FORMAT (F10.6,F10.6,F10.6) 

FORMAT ("0°,12,5X,'E =',F10.6,5X,"EK2 =',F10.6,5X, 


X'DN/DN(1) =',F10.6,5X,*EK2*DN =',F10.6) 


FORBAT (*1*,*QUPPER PLATE BIOT NO. = ',F10.6,10%, 


K'BOTTCM PLATE BIOT NO. =',F10.6) 


FORMAT {*0*,80X,810.4,5X%,£10.4) 
FORMAT (*0*,5X,*SUPERHEAT RATIO =',F10.6,10X, 


X*ICE FREE LENGTH =',F10.6) 


K FOR ICE = 1.28BTU/{HR-FT-F) 
K FOR WATER = 0.332 BTU/(HR-FT-F) 
KL/KS APPROX. 0.25 


READ (5,901) XMAX,PR 

IF (XMAX.EQ.-0.D0) GO TO 300 
BIS=BI2/4.D0 

TE=4.D0/T 

WRITE (6,888) TE,PR 
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X 1=-DELTA (1) 

PXF=1.44D2*PR*XE 

Py1=0.D0 

QT=0.D0 

T=1I+1 

IF (I-10) 5,5,6 

J=J+1 

IF {DELTA(J)) 21,21,51 

I=1 

TF 4eP;h 7,748 

D=DELTA (J) * ((1.D0/DEL1) + (1.DO0/DEL) ) /2.D0 
7p=1 

GO TO 9 

D=DELTA (J) *((1.D0/DEL2) + (4.D0/DEL1) + (1.DO/DEL)) 


X/3.0D0 


IP=0 

BN=A+D 

SUMI=0.D0 

SUMJ=0.D0 

DO 10 M=1,10 

ARG=-SQ (M) *BN 

IF (8.D14+ARG) 25,25,22 

C=EK (M) ¥DEXP (ARG) 

SUMI=SOUMI+C 

SUMJ=SUMJ- (6.D0*C) /SQ (M) 
CONTINUE 

GO TO 23 

WRITE (6,26) 

FORMAT (5X,*ARG. EXCEEDS 80.0") 
DN=SUMI* (1. D0#BIS) / (BIS* (SUMI-TE) ) 
IF (DN.LT.1.D-4) GO TO 21 

IF (DN-DFL+0.1D-3) 11,12,13 

IF (DN-DEL-0.1D-3) 12,12,11 
DEL=DN 


MAT (*4°,10X,*TW=",F10.6,10X,'PR. NO. =',F10.6) 
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K=K+1 
IF (K-15) 14,14,15 
WRITE (6,16) 
FORMAT (1X,'NOT CONVERGING") 
GO TO 21 
IF (IP) 8,8,7 
WRITE (6,1000) SUMI,SUMJ 
FORMAT (* ',30X,E14.6,5X,214.6) 
QINT=SUMI/DN 
RNU=-QINT/SOMI 
RNUU=-QINT/ (SUMJ+1.D0/T) 
IF (ISTT) 200,200,201 
INT 1=CINT 
ISTT=1 
QT=QT-3. DO*DELTA (J) * (QINT14QINT) 
DELIT=DEL1*DEL1*DEL1 
DNT=DN*DN*DN 
IF (IP) 28,28,29 
S=DELTA (J) *((1.D0/DEL1T) +(1.D0/DNT) )/2. D0 
Go TO 31 
S=DELTA (J) *((1.D0/DEL2T) + (4. DO/DEL1T) +(1.DO0/DNT) ) 


X/3.D0 


PY=1.44D2*PR* (B+S) 

X2=X1 

X1=X 

X=X4#DELTA (J) 

K=1 

DDEL= (DN-DEL2) / (X-X2) 

D2DEL= ( (DN-DEL1) / (X-X1) - (DEL1-DEL2) 7 (X1-X2) ) / (X-X2) 
PM=6.D0* (1. D0-DEL1*DEL1) / (5.D0*DEL1*DEL1) 
P=PM+PV1 

TF UP) 34,34, 27 

WRITE (6,18) X1,DEL1,DDEL,D2DEL, BN 

FORMAT (°0°,2X,F9.6,4(2X,E14.6) ) 
TB=T*(1.D0/T+SUMJ) / (1. D0+T) 

WRITE (6,32) B,PM,PV1,0T1,TB,RNU 

FORMAT (* ',4(2X,F14.5),20X,F14.5,2X,F14.5) 
XX=X1+XE 

PT=P+PXF 

COQ=QOT1+CXF 

WRITE (8,1111) X1,XX,DEL1,PT,Q00,TB,RNU,RNUD 
FORMAT ({* *',£10.3,£10.3,6 (F 10.3) ) 

DEL2=DEL1 

DEL1=DEL 

DEL=DN 

DEL2T=DEL1T 
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PV1=PV 

QINT1=CINT 

OT1=QT 

TF (TEy 020579 
A=BN 

B=B+S 

IF (X-XMAX) 4,4,21 
CONTINUE 

STOP 

END 


SUBROUTINE FF (Y,EE,C,W) 
IMPLICIT REAL*8 (A-H,0O-2) 
DIMENSION RM1(200),RM2 (200) 
E=EE/4.D0 
YY=Y-1.D0/2.D0 
RM1(1)=(1.D0-£) EEX YY*YY/2.D0 
RM2 (1) =(3.D0-E) *EE*YY*YY/6. DO 
R1=RM1(1) 
R2=RM2 (1) 
DO 10 N=2,150 
FCT=2. D0*N 
M=N-1 
RM1(N)=(4.D0*N-3.D0-E) *RM1 (M) XEEXYY*YY/ (FCT* (FCT-1.D0) ) 
RM2(N) = (4. DO*N-1.D0-E) *RM2 (M) *EEXYY*YY/ (FCT* (FCT+1. D0) ) 
R1=R1+RM1(N) 
R2=R2+RM2 (N) 
CONTINUE 
IF (DABS (RM1(150)-RM1(149)).LT.0.1D-7) GO TO 11 
IF (DABS (RM2(150)-RM2(149)).LT.0.1D-7) GO TO 11 
WRITE (6,903) 
STOP 
W= (CX (1. D04#R1) +DSQRT (EF) *YY* (1. D0+R2) ) 
X/DEXP (FE*YY*YY/2.D0) 
FORMAT ("NEED MORE THAN 150 TERMS') 
RETURN 
END 


We 


pa ; 
ona 
- 


<) g = 


as 


= 
a 


ee aa Se aes 
(W,.2 ,aa% .¥). 3 5%. aurrue asue 
(X-0,8-8) B*TAS ® TION: rai 
(O08 Say @ S)fum woreamnrg 
- . OF 0% g 
* as S\Od. P= aa 

OG .S\"Y" YY *s, Bp qr 
-O\YY* aye rs (3 Oos€) ry sa 
hin ht vase i 
or ee a be A ‘ca 
‘ORE, Sait OF Od 
“*0 S® wha 


| a. om 
; 1) \¥Y*S Yaa" (ny Phe Ch<00, E-u*oa. bys = (ML aa 
» eT vYTere nt ro) . (t-O0.f~H*0d, By = r (ul Sea 


tl ae 

(Cer yPRA- (OPTS EMAPARAG) az 

(VET) SRA (Oct) SMR pedagy ar 

(£00, 3), aPTAW 

; iF ne gore 

((S2+00 .0)*YY* (7a) THOZG (f 4ONS) 8D =e OF 
(OG. S\YY*¥ ens ay ak 

(‘SUT OF KART SHON CTEM") TARROE 9 ERE 


cc 
cc 
Se 
cc 
cc 
cc 
cc 
cc 


cc 
cc 
i 


cc 
cc 
cc 


cc 


Cc 
cc 


100 
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PROGRAMS FOR SOLIDIFICATION FREE ZONE 
INCLUDING AXIAL HEAT CONDUCTION 


1. EIGENVALUES 


IMPLICIT REAL*8 (A-H,O-Z) 
DIMENSION F (3) 

DC 100 IT=%, 12 

READ (5,901) BI,PE 
N1=150 

WRITE (6,1001) BI,PE 


EIGENVALUES 
DG 100 I=1,20 


READ (5,902) XL 
READ (5,902) XR 


FALSE POSITION 


30 


CALL HYPFN (PE,XR,N1,FR,F,8,N2) 
CALL HYPEFN (PE,XL,N1,FL,F,%,N2) 


ITERATICNS 


DO 20 KK=1,30 

XM=XL- (XR-XL) *FL/ (FR-FL) 

CALL HYPFN (PE,XM,N1,FM,F,R8,N2) 
WRITE (6,1004) KK,FM,XM 

IF (DAES (FM).LT.1.D-8) GO TO 21 
IF (FL¥FM) 1,21,3 

XR=XM 

FR=Fa 

GO TO 20 

XL=XM 

FL=Fil 

CONTINUE 

WRITE (6,1005) I,XM,N2 

WRITE (6,1002) R 

WRITE (7,903) XM,R 

CONTINUE 

FORMAT (D15.8,D15.8) 

FORMAT (L8.4) 
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FORMAT (D20.10,D20. 10) 


FORMAT ('1',10X,'BICT NO. =',F10.5,15X, 


X*PECLET NO. =',F10.5) 


FORMAT (r0", 85X,'R(1) =',D15.8) 


FORMAT (* ',12,5X, "FU =" p15.8, 10%, 
FORMAT (*0',80X,*EIGENVALUE ',12,? 
X10%) t82%=' ,73) 


STOP 
END 


2» EIGENCONSTANTS 


IMPLICIT REAL*8 (A-H,0-2Z) 
REAL*4 EPS 


—_ 


ig 
7 ¥ 


F14 


ae 8) 


DIMENSICN AL(20) ,BE(20) ,¥N(20,201) ,8N(20,201) 


DIMENSION A(40,40),B(40) ,W1(201),W2 (201) 


DIMENSION W(201),2(201),D(20,20) ,AA (1600) 


READ (5,902) (AL(I) ,I=1,20) 


READ (5,903) ((YN(N,K),K=1,201) ,N=1 


SG JO 2716.) | INFINITY 


DO 950 Ic=1,4 

READ (5,301) BI,PE 

LF UBISLT. 420170}. CO’ 40) 1 
WRITE (6,1009) PE 

GO TO 2 

WRITE {6,1001) BI,PE 

WRITE (6, 1002) 

WRITE (6,1003) ({AL(1) ,I=1, 20) 
N1=150 

DO 10 J=1,20 

READ (5,904) BE(J) 

E=BE (J) 

A=1.D0+ (6.4D1*E*E) / (9. DO*¥PE*PE) 
Al= (1. D0-A¥*E) /4.D0 

Y=-0.5D-2 

DO 10 K=1,201 

Y=Y+0.5D-2 

CALL FURY (22,E,41,%,81,88, 82) 
RN(J,K)=RR 

CONTINUE 

WRITE (6,1004) 

WRITE (6,1003) (BE(I) ,I=1,20) 
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NDIM=201 

H=0.5D-2 

DO 200 L=1,20 

DO 200 M=1,20 

DO 100 N=1,201 
W(N)=RN (M,N) *YN(L,N) 
CONTINGUE 

CALL DQSF (H,W,2,NDIM) 
A(L,™)=Z (201) 

LL=20+4L 

MM=204+M 
A(NM,LL)=AL (L) *AL (L) *A (L,®) 
CONTINUE 

WRITE (6,1011) 

WRITE (6,1008) ({A(I,J) ,J=1,20) ,I=1,20) 
WRITE (6,1012) 

WRITE (6,17008) ({A(1I,J),J=21,40) ,1T=21,40) 
CALL SYMX (H,NDIM,YN,D) 

DO 300 I=1,20 

JIJ=20 

DO 300 J=1,20 

JJ=JJ+1 

A(I,JJ)=-D (I,J) 

CONTINGE 

WRITE (6, 1012) 

WRITE (6,1008) ((A(I,d3),3=21,40) ,I=1,20) 
CALL SYMX (H,NDIAS,RBN,D) 

DO 400 J=1,20 

Tr=20 

DO 400 {=1,20 

II=1II+1 
A(T1I,J)=D (1,3) *BE (J) *BE (J) 
CONTINUE 

WRITE (6,1012) 

WRITE (6,1008) ({A(I,J),J=1,20) ,1=21,40) 
DO 600 I=1,20 

DO 500 J=1,201 

W(J)=Y¥N (I,J) 

CONTINUE 

CALL DQSF (H,W,Z,NDIM) 
B(I)=Z (201) 

CONTINUE 

DO 700 K=21,40 

B(K)=0.D0 

CONTINUE 

WRITE (6,1013) 
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WRITE (6,1008) (B(I),I=1,40) 
DO 1750 J=1,40 
DO 750 I=1,40 
IJ= (J-1) ¥404I 
AA(IJ)=A (I,J) 
750 CONTINUE 


M=40 
N=1 
EPS=1.5-14 
os 
CC SOLVE SYSTEM OF EQUATIONS 40*40 
oe A*X=E 
Ce 


CALL DGELG (B,A,M,N,EPS,IER) 
WRITE (6,1005) IER 
WRITE (6,1006) (B(I),I=1,40) 
WRITE (7,903) (B(1),1I=1,40) 
DO 900 T=1,201,190 
T1=1.D0 
T2=0.D0 
DC 800 J=1,20 
JJ=204+3 
TI=T1I+E (JJ) *YN (J,1) 
T2=T2+B (J) *RN(J,1) 
WRITE (6,1010) 71,72 
800 CONTINUE 
WRITE (6,1007) I,T1,T2 
WRITE (7,1015) T2 
900 CONTINUE 
CC 
CC CHECK BOLK TEMP. AT X=0 
ce 
Y=0.D0 
DO 920 IT=1,201 
$1=0.D0 
S$2=0.D0 
JJ=20 
DO 910 J=1,20 
JJ=JJ+1 
$1=S1+B2 (JJ) *YN (J,T) 
$2=S2+E (J) *RN (J,T) 
910 CONTINUE 
W1(I)=(1.D0-Y¥*Y) *31 
W2 (I)=(1.D0-Y¥*Y) *S2 
Y=Y+H 
920 CONTINGE 
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CALL DCSF (H,W1,Z,NDIM) 
T1IM=1.D0+1.5D0*Z (201) 
CALL DQSF (H,W2,2Z,NDIM) 
T2M=0.15D1¥*Z (201) 

WRITE (6,1014) T1M, 72M 
CONTINUE 

FORMAT (D15.8,D15.8) 
FORMAT (D16.9) 

FORMAT (5D16.9) 

FORMAT (D20.10) 

FORMAT (*1%,20X,'BIOT NO. =',F12.4,15X,"PECLET NO. =' 


X,F12.4) 


FORMAT ('0','ALPHA :! 

FORMAT gh elapse, 

FORMAT ('0',*BETA :*) 

FORMAT ('0',*EIGENCONSTANTS :',50X,'IER =',12) 
FORMAT ('0°,10X,5D16.9) 

FORMAT ('0',20X,13,5X,'T1 =',D16.9,10X,'T2 =",D16.9) 
FORMAT (' ',5D16.9) 

FORMAT (*1°,15X,'BIOT NO.= INFINITY',15X,'PECLET NO. 


X=',F12.4) 


FORMAT (' *',5X,D16.9,5X,D16. 9) 

FORMAT (*0'," : MATRIX A: °') 

BORBAT (4PO855E 8) 20 8¢ 20RET Pe +, 70K," 2 84 70%,% = 4) 
FORMAT ('0'," : MATRIX B=: ‘) 

FORMAT (10',5X,*X=0',19X,"T1IM =",D16.9,9X,'T2H =! 


~~ ue =a 6 


X,D16.9) 


FORMAT (F7.4) 
STOP 
END 


3. BULK TEMP. AND LOCAL NUSSELT NO. 
DCWNSTREAM AND UPSTREAM 


IMPLICIT REAL*8 (A-H,O-Z) 

DIMENSION AL(20) ,BE(20) ,YN(20,201) , RN (20,201) ,B (20) 
DIMENSION C(20),X(35),W(101),Z (101) ,D(101) , DR (20) 
H=0.1D-1 

NDIM=101 

N1=150 

READ (5,901) (X(I),I=1,35) 

READ (5,902) BI,PE 

WRITE (6,1001) BI,PE 

READ (5,905) (AL(I) ,1=1,20) 
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READ (5,903) (BE(I) ,I=1,20) 

READ (5,904) (C(I) ,I=1,20) 

READ (5,904) (B(I),I=1,20) 

READ (5,909) ({¥N (I,J) ,d=1,201) ,I=1,20) 
DO 200 I=1,20 


Y=-H 
DO 200 J=1,NDIM 
Y=Y+H 


A=1.D0+ (6.4D1*BE (I) *BE (I) ) / (9. DO*PE*PE) 
A1=(1.DO-A¥*BE(I))/4.D0 

CALL FUNY (PE,BE(I),A1,Y,N1,RR, N2) 
RN(I,J)=RR 

CONTINUE 

Y=1.D0 

DO 150 I=1,20 

A=1.D0+ (6.4D1*BE (I) *BE (I) ) / (9. DOXPE*PE) 
Al=(1.DO0-A*BE(I))/4.D0 

A2=A1+1.50 

CALL FUNY (PE,BE(I),A1,Y¥,N1,F1,N2) 
CALL FUNY (PE,BE(I) ,A2,Y,N1,F2,N2) 
DR (I)={ (-BE (1) -2. DO*A1) *F14#2.D0*A1*F2) 
CONTINUE 

WRITE (6,906) 

WRITE (6,1005) (AL(I),I=1,20) 

WRITE (6,1005) (B(I),I=1,20) 

WRITE (6,906) 

WRITE (6,1005) (BE(I),I=1,20) 

WRITE (6,1005) (C(I) ,I=1,20) 

WRITE (6,1004) YN(10,201) ,YN (20,201) 
WRITE (6,1007) RN(10,101) ,8N (20,101) 
WRITE (6,906) 

WRITE (6,1005) (DR(I) ,I=1,20) 

WRITE (6,906) 

DO 100 K=1,10 

Y=-H 

IF (X(K)-LT.0.D0) GO TO 51 

DO 90 KK=1,101 

Y=Y+H 

SUM1=0.D0 

DO 80 L=1,20 

S=BE(L) *BE (L) 

Q1=8. DO*S*X (K) /3.D0 

IF (Q1.GT.8.D1) GO TO 81 
SUM1=SUM14C (L) *RN (L,KK) /DEXP (Q1) 
CONTINUE 

D (KK) =S0M1 
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W (KK) =SUM1* (1. DO-¥*Y) 
90 CONTINUE 
CALL DCSF (H,W,Z,NDIM) 
TM=1.5D0*Z (101) 
SUM2=0.D0 
po 95 I=1,20 
U=BE (I) *BE (I) 
Q2=8.D0*U*X (K) /3.D0 
IF (02.61.7.D1) GO TO 96 
SUM2=SUM2+C (I) *DR (I) /DEXP (02) 
95 CONTINUE 
96 RNU=-2.D0*SUM2/TM 
WRITE (6,1002) X(K),TM,RNOU 
WRITE (6,1006) (D(I) ,I=1,101,10) 


GO TO 100 
is 
CC UPSTREAM SCLUTION 
os 
51 DO 40 #=1,201,2 
Y=Y+H 
SUM3=1.D0 


DO 30 N=1,20 
Q3=8. DO*AL (N) *AL (N) *X (K) /3.D0 
Pr 102.6787 204) GO TA37 
SUM3=SUM3+B (N) *YN (N,M) *DEXP (C3) 
30 CONTINUE 
31 NN= (M41) /2 
W (NN) =SUM3* (1.D0-Y*Y) 
D (NN) =SU63 
40 CONTINUE 
CALL DQSF (H,W,Z,NDIM) 
TM=1.5D0*Z (101) 
WRITE (6,1003) X(K),TM 
WRITE (6,1006) (D(1I),1I=1,101,10) 
100 CONTINUE 
901 FORMAT (5D10.4) 
902 FORMAT (D10.4,D10.4) 
903 FORMAT (020.10) 
904 FORMAT (5D16.9) 
905 FORMAT (D16.9) 
906 FCRMAT (*'0',/) 
1001 FORMAT ("1*,15%,"BIOT NO. =',F10.4,10X,"PECLET NO. 
X,F10.4) 
1002 FORMAT (°0',5X,'X ="',F8.4,5X,'*BOULK TEMP. =*',F10.4, 
X10X,"NU NO. =',F15.4) 
1003 FORMAT ('0',5X,*X =',¥F8.4,5X,"BULK TEMP. =',F10.4) 
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1004 FORMAT ('0*,5X,"YN (10,201) =",D16.9,10X, "YN (20,201) =! 
X,D16.3) 
1005 FORMAT (5X,5D16.9) 
1006 FORMAT (* *,80X,F10.4) 
1007 FORMAT (*0*,5X,'RN(10,101)=',D16.9,10X, RN (20,101)=' 
X,D16.9) 
STOP 
END 
cc 
i, 
SUBROUTINE HYPPN (PE,E,N1,F0,F ,8,N2) 
IMPLICIT REAL*8 (A-H,O-Z) 
DIMENSION F (3) 
A=1.D0+4 (6.4D1*E*E) / (9. DO*PE*PE) 
Al=(1.00-A*E) s4.D0 
NAI=A1 
IF (A1.GE.-2.D0) GO TO 14 
cc 
ec xe F(1)= M (AJ*+2 , 1/2 , EB ) *** 
cc Kee F(2)= M (AJ+1 , 1/2 , E ) *** 
i 
TI=A1-NA1 
T2=T1+1.D0 
CALL FCN (PE,T2,E,N1,F(1) ,DUMP,N2) 
CALL FCN (PE,T1,E,N1,F(2) ,DUMP,N2) 
cc 
te 4*%* GENERATE F(3)= M (AJ, 1/2 , E) *** 
cc 
M1=-NA1 
DO 50 I=1,81 
T1=T1-1.D0 
F (3) =- ((114+1.D0) *F (1) - (2.D0*T14#E+1.5D0) *F (2) 
X/ (T14+0.5D0) 
F (1) =F (2) 
F (2) =F (3) 
50 CONTINUE 
FU= (BI-E-2.D0*A1) *F (2) +2. D0*A1*F (1) 
R=F (2) /DEXP (E/2. D0) 
GO T0 61 
44 CALL FCN (PE,A1,E,N1,FN,FU,N2) 
R=FN/DEXP (E/2.D0) 
DO 60 J=1,3 
F (J)=0.D0 
60 CONTINUE 
61 RETURN 
END 
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SUBROUTINE FCN (PE,T,E,N1,FN,FU,N2) 


*HRHHEN T.LT.-2 THIS SUBROUTINE 
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IMPLICIT REAL*8 (A-H,O-Z) 
DIMENSION RM1(400) 
FU=0.D0 

N2=N1 

AM=4.D0*T 

RN1 (1) =AN*E/2.D0 
FN=1.D0+8™1 (1) 

FCT=2.00 

DO 10 L=1,N2 
FCT=FCT#+2.D0 

AM=AM4+4.D0 

M=L41 

RM1(@%) =AM*RM1(L) *E/ (FCT* (FCT-1. D0) ) 
FN=FN+R81(M) 

CONTINUE 


IF ( (DABS ({RH7(N2+4) ) +DABS (RH 1(N2))).LT.7.D-12 ) 


XGC TO 12 

N2=N2450 

IF (N2.GE.400) GO TO 11 
GC TO 9 

WRITE (6,1003) 

STOP 

IF {T.L7.-2.D0) GO TO 37 
FU=BI-E 

DO 30 K=1,N2 
FU=FU4+RM1(K) * (BI-~E+2.D0*XK) 
CONTINUE 

FORMAT ("NEED MORE THAN 400 TERMS') 
RETURN 

END 


SUBROUTINE FUNY (PE,E,A1,Y,N1,R,N2) 
IMPLICIT REAL*8 (A-H,O-Z) 

DIMENSION F (3) 

NAI=A1 

IF (A1.GE.-2.D0) GO TO 14 
T1=A1-NA1 

T2=T14+1.D0 

CALL FNY (PE,T2,E,Y,N1,F (1) ,N2) 
CALL FNY (PE,T1,£,Y,N1,F (2) ,N2) 
M1=-NA1 
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DO 50 I=1,M1 
T1=T1-1.D0 
F (3) =- ( (7141.D0) *F (1) - (2. DO¥T14+E8*Y *¥4+1.5D0) *F (2)) 
KX/ (T14+0.5D0) 
F (1) =F (2) 
F (2)=F (3) 
50 CONTINUE 
R=F (2) /DEXP (EXY*Y/2. D0) 
GC TO 61 
14 CALL FNY (PE,A1,£,Y¥,N1,FN,N2) 
R=FN/DEXE (E*Y*Y/2. D0) 
61 RETURN 
END 


SUBROUTINE FNY (PE,T,E,Y,N1,FN,N2) 
IMPLICIT REAL*8 (A-H,O-Z) 
DIMENSICN RM1(400) 
N2=N1 
9 AM=4.D0*T 
RM1(1) =AM*EXY*Y/2.D0 
FN=1.D04+RM1(1) 
FCT=2.50 
pO 10 L=1,N2 
FCT=FCT+2.D0 
AM=AM44.D0 
M=L+1 
RM1 (M) =AM*RM1(L) ¥E*Y*Y/ (FCT* (FCT-1. D0) ) 
FN=FN4RH1(M) 
10 CONTINUE 
IF ( (DABS (RM1(N2+41))+DABS (RM1(N2))).LT. 1.D-12) 
XGO TO 12 
N2=N2+50 
IF (N2.GE.400) GO TO 11 
GO TO 9 
11 WRITE (6,1003) 
STOP 
12 CONTINUE 
1003 FCRMAT ("NEED MORE THAN 400 TERMS*) 
RETURN 
END 
ee 
SUBROUTINE SYMX (H,NDIM,FN,D) 
IMPLICIT REAL*8 (A-H,O-2Z) 
DIMENSION FN(20,201),D(20,20),W(201),Z(201) 
K=0 
pO 2 I=1,20 
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K=K+1 

DO 2 J=1,K 

DO 1 N=1,201 
W(N)=FN(I,N) ¥FN(J,N) 
CONTINUE 

CALL DCSF (H,W,2Z,NDIM) 
D (I,J) =2Z (201) 

IF (1.80.3) 60 TO 2 
D(d,T)=D) (1,3) 
CONTINUE 

RETURN 

END 
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